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1  Objective 

The  goal  of  this  project  was  to  expand  existing  investigations  and  processing  routines  for  so- 
called  “Advanced  EMI”  arrays.  Previous  demonstrations  conducted  with  advanced  systems  both 
at  seeded  sites  (e.g.,  APG  and  YPG)  and  live  sites  such  as  Camp  Sibert,  AL,  Camp  San  Luis 
Obispo  (SLO),  and  Camp  Butner  have  shown  that  unanswered  questions  remain  relating  to 
antenna  array  configuration,  physics-based  target  modeling,  and  classification.  Our  objective 
here  has  not  been  to  develop  new  sensors,  sensor  hardware,  or  programs  for  target  modeling. 
Rather,  it  is  to  seize  an  opportunity  to  review  and  update  earlier  research  in  retrospect  in  order  to 
search  for  a  better  and/or  simpler  antenna  configuration,  and  to  maximize  the  benefit  derived 
from  these  systems  through  expanded  data  processing  and  better  anomaly  classification.  The 
study  we  conducted  was  enabled  and  justified  by  the  fact  that  real  data  and  ground  truth  was 
available  from  the  ESTCP  classification  studies  at  SLO  and  Camp  Butner  in  North  Carolina. 
The  research  has  been  focused  on  five  specific  objectives: 

1.  Improved  Target  Detection  -  The  objective  of  this  task  was  to  exploit  the  capability  of 
an  instrument  such  as  the  MetalMapper,  which  samples  the  vector  field  at  multiple 
locations  inside  the  transmitter  loop,  to  generate  simple  physics-based  estimates  of 
position  and  size  of  targets  within  the  field  of  view  of  the  transmitter.  Applied  to 
dynamic  data  either  in  real-time  during  acquisition  or  during  post-acquisition 
processing,  this  type  of  processing  can  serve  as  a  primary  target  picking  screen  that 
sorts  targets  into  at  least  2  and  perhaps  3  categories:  1)  high-confidence  clutter;  2) 
targets  selected  for  cued  ID;  and  perhaps  3)  high-confidence  targets  of  interest  (i.e. 
dig).  The  advantage  here  is  to  reduce  the  number  of  anomalies  that  must  be  revisited 
for  static  data  collection.  These  algorithms  have  been  implemented  as  part  of  the 
acquisition  software,  allowing  the  operator  to  see  a  higher-level  processed  answer  in 
real-time  to  assist  in  re-locating  static  targets,  or  to  stop  and  acquire  static  data. 

2.  Target  Parameter  Extraction  with  Multiple-Point  Static  or  Dynamic  Data  Sets  -  The 
MetalMapper  and  other  advanced  EMI  systems  were  developed  with  cued  target 
characterization  in  mind.  In  principle,  a  data  set  acquired  at  a  single  spatial 
measurement  point  is  sufficient  for  the  target  analysis.  Here,  our  objective  has  been  to 
determine  the  extent  to  which  the  quality  of  parameter  extraction  is  enhanced  by 
sampling  the  secondary  EMI  response  at  several  spatial  points  in  the  proximity  of  the 
target.  In  addition,  we  have  analyzed  the  ability  to  discriminate  targets  using  the 
dynamic  data.  We  find  that  some  of  the  benefits  obtained  from  more  sample  points 
over  a  larger  aperture  will  be  offset  by  two  mechanisms.  First,  the  platform  is  in 
motion,  producing  noise  and  position  uncertainty.  Second,  the  dynamic  data  acquired 
with  the  MetalMapper  used  only  the  vertical  axis  transmitter. 

3.  Optimal  Array  Configuration(s)  -  The  goal  was  to  define  one  or  more  antenna  arrays 
that  meet  the  need  for  cost-effective  sensor  arrays  for  static  (i.e.  cued)  measurements 
for  characterization  and  identification.  This  could  result  in  smaller  systems  that  are 
easier  to  deploy  and  less  expensive,  yet  may  provide  nearly  equivalent  performance 
for  discrimination. 
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4.  Multiple  Source  Detection  -  The  objective  here  is  to  identify  those  target  picks 
associated  with  complex  secondary  magnetic  fields  that  arise  either  from  the  presence 
of  multiple  targets  in  close  proximity,  targets  that  are  asymmetrical,  or  that  may  be  too 
large  dimensionally  to  be  approximated  with  a  point  dipole.  A  multi-target  inversion 
algorithm  has  been  developed  and  tested.  Downward  continuation  and  reduction-to- 
pole  algorithms  have  also  been  developed  to  aid  interpretation  and  to  provide  useful 
starting  models  for  the  multi-target  inversion. 

5.  Oasis  montaj  Integration  -  This  task  was  not  part  of  the  original  work  scope,  but  was 
added  after  the  program  office  suggested  that  the  new  routines  should  be  integrated 
into  existing  Oasis  montaj  workflows. 

Although  some  of  these  objectives  overlap  (e.g.  detecting  multiple  targets),  the  methods  for 
reaching  these  objectives  are  different. 

Work  was  performed  on  this  project  by  Geometries  ,  Snyder  Geoscience  (SGI),  Earth  Science 
Systems  (ESS),  and  G&G  Sciences.  SGI  was  primarily  involved  with  supplying  data,  creating 
data  format  conversion  routines,  testing  the  codes,  and  conducting  the  optimal  array 
configuration  analysis.  ESS  wrote  the  software  for  modeling  data,  inverting  data,  target 
detection,  complex  source  treatment,  and  Oasis  montage  integration.  SGI  and  ESS  charged  time 
against  Task  2,  “Parameter  Extraction  Using  Multi-Point  Data  Sets”  and  towards  Task  3, 
“Optimal  Array  Configurations”.  ESS  charged  time  against  Task  2,  “Parameter  Extraction  Using 
Multi-Point  Data  Sets”,  Task  4,  “Multiple  Source  Detection”,  and  Task  5,  “ Oasis  montaj 
Integration”.  Task  1,  “Improved  Target  Detection”,  was  tackled  by  G  &  G  Geosciences. 
Geometries  performed  the  overall  management  of  the  project. 
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2  Background 

With  ESTCP  support  (MM-0603),  Geometries  undertook  the  commercialization  of  an  advanced 
EMI  system,  called  the  MetalMapper,  for  application  to  UXO  detection  and  characterization.  In 
2009,  the  MetalMapper  participated  in  a  live-site  demonstration  at  the  former  Camp  San  Luis 
Obispo  (SLO)  in  central  California  (Figure  2.1),  and  in  2010  at  the  former  Camp  Butner  in  North 
Carolina.  The  data  acquired  at  SLO  and  Camp  Butner,  together  with  the  ground-truth  that 
resulted  from  the  digging  provides  a  unique  and  timely  opportunity  for  retrospective  analyses 
directed  toward  improving  the  performance  of  this  new  generation  of  EMI  sensors 


This  study  was  performed  as  a  result  of  a  desire  to  improve  data  analysis  methods  using  the 
lessons  learned  during  the  MetalMapper  surveys  at  San  Luis  Obispo  (SLO)  and  Camp  Butner. 
ROC  curves  from  those  surveys  (Figure  2.2  and  Figure  2.3)  are  excellent,  but  we  wish  to 
investigate  the  failures  to  further  understand  the  limitations  and  improve  the  performance  for  the 
small  portion  of  UXO  that  are  labeled  as  non  TOI. 


3 


Expanded  Processing  Techniques  for  EMI  Systems 


SERDP  MR- 1772 


ROC  Cuives  (Test):  SLO  Discrimination 
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Figure  2.2:  Performance  curves  from  San  Luis  Obispo 
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Figure  2.3:  Performance  curves  from  Camp  Butner 
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3.1  MMTargets 


3  Materials  and  Methods 


A  new  modeling  and  inversion  program  called  MMTargets  was  created  for  this  project.  Derived 
from  an  existing  program  that  was  written  for  the  USGS  ALLTEM  system  to  invert  dynamic  data 
[1],  MMTargets  can  be  used  to  generate  synthetic  datasets,  or  invert  data  to  estimate  target 
parameters.  The  program  has  a  large  number  of  user-configurable  settings  that  control  how 
synthetic  data  are  generated,  or  how  the  program  processes  and  inverts  a  dataset.  These  settings 
are  configured  in  an  INI  file. 

MMTargets  can  model  both  static  and  dynamic  data.  In  the  forward  modeling  mode,  it  provides 
options  to  add  simulated  EM  noise  to  the  synthetic  data.  Different  noise  levels  can  be  added  for 
vertical  and  horizontally  oriented  coils.  Uncertainty  in  positional  and  orientation  are  simulated 
by  adding  random  noise  to  those  parameters  before  generating  data.  The  inversion  routine  uses 
several  constrained  models  when  inverting  data.  The  inverted  datasets  can  be  a  single  cued 
snapshot,  multiple  queued  snapshots,  or  dynamically  acquired  data.  Furthermore,  several 
different  target  shape  constraints  can  be  enforced;  and  the  results  provide  better  information  for 
classification  routines. 

MMTargets  provides  flexible  and  highly  configurable  inversion  routines.  In  addition  to  target 
shape  constraints,  limits  can  be  placed  on  the  location  of  the  target  and  its  maximum  size.  These 
limits  are  typically  set  using  the  extents  of  the  survey  coverage  and  the  size  limits  of  the  targets 
at  a  given  site.  The  target  location  can  be  constrained  using  the  size  and  location  (centroid)  of 
the  anomaly,  or  through  anomaly  focusing  efforts  such  as  reduction-to-pole  or  downward 
continuation  (discussed  below).  The  general  inversion  sequence  is  as  follows: 

1.  Determine  a  reasonable  set  of  initial  parameters  by  inverting  a  simple  4-parameter 
spherical  model  to  estimate  location  and  sphere  diameter. 

2.  Perform  a  constrained  inversion  using  any  number  (specified  by  user)  of  the 
following  models:  sphere,  oblate  spheroid,  prolate  spheroid,  ellipsoid.  Users  can  also 
specify  bounding  constraints  for  each  model  parameter.  This  is  useful,  for  instance, 
when  the  lateral  location  of  the  target  is  well  known. 

3.  Select  new  starting  models  from  a  list  or  optionally  generate  new  starting  models 
(Monte  Carlo,  genetic  mutation)  and  repeat  the  constrained  inversion. 

4.  Generate  principal  polarizability  curves  by  separately  diagonalizing  the  data  at  each 
time  gate.  Insure  that  each  polarizability  curve  has  directional  continuity  in  time  so 
that  principal  polarizability  curves  do  not  switch  at  cross-overs. 

The  principal  polarizability  curves  produced  by  MMTargets  are  not  a  result  of  the  joint 
diagonalization  procedure  that  is  commonly  applied  in  UXO  processing.  Because  the  Eigen 
currents  precess  as  the  currents  decay,  estimating  a  rotation  matrix  that  best  diagonalizes 
polarization  matrixes  for  all  times  is  not  a  fair  representation  of  the  physical  process.  An  average 
Eigen-structure  does  not  account  for  precession.  When  MMTargets  generates  polarizability 
curves  it  diagonalizes  the  polarizability  matrix  at  each  time-gate,  and  insures  that  there  is 
directional  continuity  in  each  curve  from  time-gate  to  time-gate.  The  usefulness  of  this  method 
is  illustrated  by  re-processing  the  data  from  targets  at  the  Camp  Butner  demonstration  that  were 
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missed  in  earlier  processing.  Target  849  (from  our  anomaly  list,  program  office  number  404)  is  a 
very  typical  37  mm  ordnance,  but  the  polarizability  curves  generated  by  the  inversion  program 
MMRMP  suggested  a  plate-like  object,  while  those  from  MMTargets  suggest  a  rod-like  object 
(see  Figure  3.1). 


MMTargets  MMRMP 

Figure  3.1:  Comparison  of  the  polarizability  curves  produced  by  MMTargets  and  by  MMRMP.  MMTargets  results 
correctly  suggest  a  rod-like  object,  while  the  MMRMP  results  suggest  a  plate-like  object 


3.2  Constrained  Single-Target  Inversion 

MMTargets  has  the  capability  to  employ  constraints  on  its  search  to  find  a  solution  that  best  fits 
the  data.  The  user  can  set  constraints  that  will  restrict  the  resulting  target  solution  to  be  a  sphere 
(SP),  a  prolate  spheroid  (PS),  an  oblate  spheroid  (OS),  or  an  ellipsoid  (EL).  The  ellipsoidal 
solution  is,  of  course,  the  solution  to  the  inverse  modeling  that  has  no  constraints.  How  can  we 
use  these  constraints  as  an  aid  in  discrimination?  In  this  section,  we  will  suggest  how  the  results 
from  a  constrained  solution  can  be  used  to  form  an  F-Test  statistic  the  can  be  used  to  make  a 
decision  about  the  symmetry  properties  of  a  target. 

In  Section  4.3  we  show  an  example  (see  Figure  4.9)  of  how  the  classic  expression  of  symmetry 
for  a  buried  target  (i.e.,  a  major  polarizability  curve  and  two  smaller  but  nearly  identical  minor 
polarizability  curves)  can  be  compromised  simply  by  the  fact  that  the  MetalMapper  platform  is 
laterally  offset  from  the  target.  In  the  case  shown  in  Figure  4.9,  we  were  lucky  to  have  taken  a 
repeat  measurement  that  was  located  more  directly  over  the  target  and  the  resulting  extraction  of 
the  principal  polarizability  curves  exhibited  the  classic  behavior  of  an  elongated  body  of 
revolution.  If  we  don’t  have  a  repeat  measurement  that  is  well  centered  over  the  target,  can  we 
develop  a  feature  that  is  more  robust  than  the  symmetry  feature  (Poy/Poz)  we  used  in  the  optimal 
array  study?  In  this  report  we  show  that  the  answer  to  this  question  is  YES. 

The  MMTargets  single-target  inversion  provides  a  rich  set  of  properties  that  can  be  supplied  to 
classifiers.  For  each  model  (sphere,  prolate  spheroid,  oblate  spheroid,  and  ellipsoid)  the 
estimated  parameters  are  the  target  location,  principal  radii,  orientation,  and  the  data  misfit.  All 
of  these  parameters  provide  valuable  information  for  classification  routines.  For  example,  if  the 
data  misfit  for  the  spherical  model  and  the  prolate  model  are  similar,  and  the  prolate  model  has  a 
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large  aspect  ratio,  this  reduces  the  confidence  that  the  target  actually  has  a  large  aspect  ratio.  If 
multiple  models  provide  similar  misfit  metrics,  this  may  indicate  that  a  complex  (non-dipole) 
target  is  present  and  more  advanced  analysis  algorithms  may  be  needed. 

3.3  The  F-Test  Statistic 

Using  the  capability  of  MMTargets  to  fit  the  observed  data  with  4  different  models,  3  of  which 
are  constrained  to  be  bodies  of  revolution  (BOR)  we  can  form  a  statistic  similar  to  an  F-Test  as  it 
is  traditionally  applied  to  regression  problems.  The  statistic  is  a  relationship  between  the  mean 
square  errors  of  the  constrained  solution  (MSEc)  and  that  of  the  unconstrained  solution  (MSEel)- 
For  the  purposes  of  this  section  we  shall  assume  that  MSEel  indicates  the  MMTargets  solution 
for  the  ellipsoidal  model.  We  will  use  MSEsp,  MSEps,  and  MSEos,  respectively,  according  to 
whether  the  constrained  solution  corresponds  to  the  sphere  (SP),  prolate  spheroid  (PS),  or  the 
oblate  spheroid  (OS).  We  have  experimented  with  various  relationships  for  the  F-Test  and  the 
one  that  works  the  best  for  our  purposes  is 

p  _  MSEc  —  MSE  EL 
c  MSEel 

Presumably,  the  unconstrained  solution  will  always  have  an  MSE  that  is  smaller  than  those 
generated  by  a  constrained  solution.  Thus  Fc  is  always1  a  number  that  is  greater  than  0.  In 
words,  the  statistic  is  simply  a  measure  of  how  much  larger  the  constrained  solution  error  is 
relative  to  the  fully  unconstrained  solution.  Basically,  the  idea  is  that  when  the  Fc  statistic  for 
the  constrained  solution  is  close  to  zero,  it  suggests  that  we  cannot  discount  the  possibility  that 
the  target  has  the  shape  indicated  by  the  constraint  in  question.  Thus,  the  Fc  statistic  can  be 
viewed  as  an  indicator  of  target  symmetry  characteristics. 

3.4  F-Test  -  Synthetic  Target  Populations 

Using  the  synthetic  data  set  for  the  60mm  prolate  spheroid  target  (i.e.  250  targets  randomly 
distributed  below  the  MetalMapper  over  an  offset  of  0.5  meters,  see  Section  4.3),  we  have 
computed  Fps,  Fos,  and  Fsp  according  to  Equation  3.1.  Figure  3.2  shows  the  distribution  of  these 
three  F-Test  statistics  for  spatially  identical  distributions  of  a  180mm  oblate  spheroid  (1:3  aspect 
ratio,  Figure  3.2A)  and  a  60mm  prolate  spheroid  (3:1  aspect  ratio,  Figure  3.2B).  The  row 
parameter  is  the  model  constraint  moving  from  prolate  spheroid  model  on  the  top  to  sphere  on 
the  bottom  row.  This  distribution  easily  shows  the  identity  of  the  target  population  as  a  whole 
since  only  the  constraint  that  corresponds  to  the  model  type  (i.e.,  oblate  spheroid/prolate 
spheroid)  produces  any  significant  numbers  of  targets  from  within  the  population  (250  targets) 
that  have  F-Test  <  0.1  (LogioFTest  <  -1).  This  is  a  useful  observation  since  it  provides  a  useful 
threshold  on  values  of  the  F-Test  statistic  below  which  we  may  find  the  F-Test  to  be  a  useful  and 
robust  indicator  of  body  of  revolution  symmetry. 


1  In  practice,  we  occasionally  find  a  situation  where  MSEc<MSEel. 
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In  Figure  3.3  below,  we  show  a 
number  of  Log-Log  scatter  plots  in 
which  the  F-Test  statistic  for  the 
prolate  spheroid  constraint  has 
been  plotted  as  the  abscissa.2  The 
ordinates  are  other  useful  target 
attributes.  With  regard  to  the  F- 
Test,  points  plotting  at  values  of 
Logio(FTest)<-l  suggest  that  the 
data  for  that  particular  target  are  fit 
well  by  a  model  constrained  to  be  a 
prolate  spheroid  (elongated  along 
the  axis  of  symmetry).  We  have 
expanded  the  scale  of  the  F-Test 
(horizontal  scale)  to  emphasize  the 
scatter  of  F-Test  values.  Two 
columns  of  figures  are  shown.  In 
the  left-hand  column  are  plotted 
points  for  all  250  targets  in  the 
target  space  (see  Figure  4.1 1).  The 
plots  in  the  right-hand  column 
contain  only  points  falling  within 
the  reduced  target  space.  In 
Section  4.3,  we  noted  that  the 
quality  of  the  parameters  extracted 
from  targets  located  outside  the 
reduced  target  space  (i.e.  targets 
with  offsets  greater  than  0.4  meters,  see  Section  4.3)  is  degraded.  This  degradation  in  the  quality 
of  the  parameter  pairs  (as  indicated  by  the  point  scatter)  is  manifest  in  the  plots  here  as  well. 
Note  that  in  the  scatter  plots  in  the  top  row  of  Figure  3.3,  both  parameters  are  indicators  of  rod¬ 
like  symmetry  and  in  a  noise-free  case  all  the  points  would  be  located  at  (F-Test=0,  Poy/Poz=1)- 

In  Figure  3.4A  we  have  taken  Figure  3.3A  (i.e.,  top  row)  and  added  corresponding  plots  based  on 
the  F-Test  statistic  for  the  oblate  spheroid  (Fos),  and  the  sphere  (Fsp)  constraints.  As  one  might 
expect  after  studying  the  histograms  in  Figure  3.2,  the  plots  based  on  Fos  and  Fsp  (Figure  3.3B 
and  Figure  3.3C,  respectively)  show  that  these  2  constrained  solutions  have  fit  errors  (MSE)  that 
are  more  than  10%  (0.1)  higher  than  the  unconstrained  (ellipsoidal)  solution. 
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Figure  3.2:  Histograms  showing  the  distribution  of  F-Test  values  for  3 
constrained  solutions  generated  by  MMTargets .  The  targets  are  a  180 
mm  oblate  spheroid  (1:3  aspect  ratio),  and  a  60mm  prolate  spheroid 
(3:1  aspect  ratio).  The  constraints  are  prolate  spheroid  (PS),  oblate 
spheroid  (OS),  and  sphere  (SP). 


2  We  focus  here  on  the  F-Test  (prolate  spheroid  constraint)  because  we  are  usually  seeking  to  identify  elongated 
bodies  of  revolution  (BOR). 
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Figure  3.3:  Scatter-plots  showing  the  behavior  of  the  F-Test  statistic  (prolate  spheroid  constraint)  when  plotted  with 
other  useful  target  features  as  the  ordinate.  The  top  row  (A)  uses  a  symmetry  feature,  the  middle  row  (B)  uses  the 
time  persistence  of  the  major  principal  polarizability  curve  (dPx/dt),  and  the  bottom  (C)  is  a  measure  of  the  aspect 
ratio  of  rod-like  body.  The  left-hand  column  shows  data  plot  for  the  full  target  set  (250  targets).  The  right-hand 
column  uses  a  reduced  target  set  (193  targets)  wherein  targets  with  offsets  greater  than  0.4m  have  been  excluded. 
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Figure  3.4:  Scatter  plot  of  a  target  attribute  indicating  rod-like  symmetry  (Poy/POz)  against  the  F-Test  attribute  with 
prolate  spheroid  constraint  (A),  the  oblate  spheroid  constraint  (B),  and  the  sphere  constraint  (C).  The  test  data  set 
consists  of  250  identical  60mm  prolate  spheroids  at  random  locations  and  random  attitudes  (see  Section  4.3  of  this 
report).  The  left  hand  column  plots  the  full  target  set  (FS).  In  the  right  hand  column,  we  plot  a  reduced  target  set 
(RS)  that  generally  excludes  targets  with  horizontal  offsets  (r)  greater  than  0.4m. 
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3.5  F-Test  and  the  Detection  of  Rod-Like  Targets 

The  distributions  that  we  display  in  Figure  3.2  suggest  that  we  can  classify  a  population  of  the 
same  target  type  (e.g.,  rod-like  targets  such  as  prolate  spheroids  or  plate-like  targets  such  as 
oblate  spheroids)  by  looking  at  distribution  of  the  F-Test  for  the  appropriate  constraint. 
However,  it  would  be  more  useful  if  the  F-Test  can  be  used  to  detect  individual  rod-like  (or 
plate-like)  targets.  The  distribution  of  F-Test  values  covers  a  multi-decade  region  of 
Logio(FTest)  space,  so  the  F-Test  cannot  be  used  in  a  2-D  scatter  plot  to  help  identify  a  cluster  of 
rod-like  bodies  of  revolution  with  similar  characteristics.  However,  generally  speaking,  the  F- 
Test  for  a  particular  constraint  (e.g.,  prolate  spheroid)  is  minimum  when  compared  with  the  F- 
Tests  for  the  two  alternate  constraints  (oblate  spheroid  and  sphere)  when  the  target  has  “rod-like” 
symmetry.  In  this  section  we  use  the  F-Test  as  a  detector  of  rod-like  bodies. 
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Figure  3.5:  Figure  showing  the  distribution  of  SNR  (A)  and  the  cumulative  distribution  of  SNR  (B)  for  750  targets 
equally  distributed  between  250  60mm  prolate  spheroids  (PS-3:1  aspect),  250  180mm  oblate  spheroids  (OS-l:3 
aspect),  and  250  100mm  spheres. 


To  demonstrate  how  such  a  detector  might  be  applied,  we  use  a  population  of  3  target  types:  a) 
60mm  prolate  spheroids  (PS);  b)  180  mm  oblate  spheroids  (OS);  and  c)  100mm  spheres  (SP). 
There  are  250  of  each  of  the  target  types  distributed  in  an  octant  of  a  spheroid  centered  with  the 
MetalMapper  array  (3T7C-MM).  The  distribution  of  these  targets  is  described  in  Section  4.3 
(see  Figure  4.11).  We  added  random  noise  to  the  model  data  consistent  with  noise  levels 
observed  at  SLO  and  reduced  by  1/V90  to  account  for  stacking.3  Figure  3.5A  is  a  stacked 
histogram  of  the  SNR  for  the  entire  target  population  (750).  Figure  3.5B  shows  the  cumulative 
distribution  for  the  3  targets.  The  colors  refer  to  each  of  the  3  different  target  types.  The  figure 


3  The  standard  deviations  of  the  noise  were  cx=30pT/s,  oy=30pT/s,  and  oz=15pT/s.  This  a  factor  of  10  less  than  we 
used  for  much  of  the  model  data  in  our  array  study.  But  in  that  study  we  did  not  account  for  a  stack  count  of  90. 
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shows  that  virtually  all  the  targets  have  an  SNR  >10 
(20dB).  Therefore,  any  inability  to  distinguish  target 
shape  from  the  target  data  sets  used  in  this  exercise 
cannot  be  blamed  on  poor  SNR. 

It  is  also  interesting  to  show  the  distribution  of  the  F- 
Test  (Fp.s-prolate  spheroid)  for  the  3  target  types.  We 
have  done  that  in  Figure  3.6.  This  distribution 
confirms  what  we  saw  previously  in  Figure  3.2.  Fps 
covers  a  multi-decade  domain  in  log-space  and  the 
distribution  of  Fps  for  the  prolate  spheroid  is 
noticeably  more  dominant  for  values  of  Fps<-1  while 
the  distribution  of  values  for  the  other  two  constraints 
(OS  &  PS)  peaks.  The  difference  is  most  obvious 
between  prolate  spheroid  (PS-Blue)  and  oblate 
spheroid  (OS-Green)  targets. 

As  a  simple  detector  of  rod-like  targets,  we  performed  the  following  rule-based  test: 

1.  For  each  target 

2.  Sort  the  4  solutions  according  to  their  respective  F-Test  values.4  In  all  but  a  few  isolated 
cases,  this  will  place  the  EL  solution  at  the  top  of  the  list  of  4  F-Test  values.  Now,  the 
value  F 2c  in  the  second  position  will  have  the  minimum  value  of  the  3  non-zero  Ftc  s. 

3.  If  F2C  <  Fth  ,  then  the  target  type  is  C  else  type  is  EL  (ellipsoid  -  meaning  not  analyzable 

for  shape).  Fth  is  a  threshold  value  for  the 
F-Test  above  which  we  designate  the  target 
as  not  analyzable  and  default  to  the  EL 
solution. 

4.  Loop  to  1  until  done. 

We  applied  the  algorithm  above  to  each  of  the  750 
targets.  As  ground-truth,  we  identified  prolate 
spheroids  (PS)  as  the  positive  decision  (  1  ),  and 
oblate  spheroids  and  spheres  and  negative  decisions 
(0).  We  used  the  fit  statistic  (MSE)  value  as  means 
of  ordering  each  decision.  That  is,  a  solution  in 
which  F 2c  corresponds  to  the  PS  constraint  and  has 
a  low  value  for  the  MSE  of  the  solution  is 
considered  a  rod-like  object  with  high  confidence. 

Likewise,  if  F2C  corresponds  to  either  an  OS  or  SP 
with  a  low  value  of  the  solution  MSE  the  target  will 
be  considered  to  be  NOT  PS  with  high  confidence. 

Figure  3.7  is  a  ROC  showing  the  performance  of  the 
F-Test  for  detection  of  rod-like  objects.  We 
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Figure  3.7:  ROC  curve  showing  the 
performance  of  the  F-Test  (PS)  for  detection  of 
rod-like  bodies  of  symmetry  (BOR).  The 
parameter  is  an  F-Test  threshold  value.  All 
targets  with  F-Test  values  larger  than  the 
threshold  are  designated  as  non-detectable 
(ND).  The  offsets  of  the  2  curves  indicate  the 
number  ND’s  found  for  the  indicated  threshold. 
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Figure  3.6:  Stacked  histogram  showing  the 
distribution  of  FPS  (F-Test  PS)  for  the  750-target 
test  population  (250  PS,  250  OS,  250  SP).  The 
colors  identify  the  type. 


4  We  set  the  F-Test  for  the  unconstrained  (EL)  solution  to  0. 
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evaluated  the  ROC  for  2  different  values  of  the  F-Test  threshold.  The  results  show  that  we  can 
detect  approximately  75%  of  the  prolate  spheroid  targets  with  a  false  alarm  rate  of  P/p=0.29. 
Most  of  the  false  alarms  are  associated  with  targets  that  are  not  detectible  (ND)  because 
FTest>FfH ■  We  have  not  calculated  the  corresponding  ROC  for  the  reduced  target  set.  It  will  be 
interesting  to  see  whether  the  false  alarm  rate  is  reduced  significantly  in  that  case.  It  is  also 
worthy  of  note  that  reducing  the  F-Test  threshold  (from  1.0  to  0.1)  does  not  appreciably  change 
the  performance  of  the  basic  rod-like  prediction.  The  only  significant  difference  between  the 
two  ROC  curves  is  the  offset  due  to  the  need  to  designate  79  more  targets  as  ones  that  we  cannot 
detect. 

The  primary  reason  for  the  broad  F-Test  distributions  is  noisy  data,  which  produces  a  noisy 
objective  function.  While  the  noise  is  small  enough  to  permit  reasonable  inversion  results,  the 
noise  causes  a  number  of  small  local  minima  in  the  same  vicinity  of  the  reported  solution.  The 
local  minima  in  the  vicinity  of  the  solution  all  have  similar  MSE  values,  and  some  will  likely 
have  lower  MSE  values  than  that  of  the  reported  solution.  It  is  anticipated  that  smoothing  the 
objective  function  will  reduce  these  difficulties  and  reduce  the  breadth  of  the  F-Test 
distributions.  This  is  recommended  as  future  work. 

More  experimentation  is  required  to  perfect  the  procedures  for  exploiting  the  information  in  the 
F-Test  values.  Clearly,  we  need  to  evaluate  the  ROC  for  larger  values  of  the  F-Test  threshold. 
And  we  also  need  to  evaluate  its  use  on  real  data.  We  have  the  data  sets  (e.g.,  SLO,  and  Butner). 
But  we  need  to  invert  each  data  set  with  MMTargets  before  we  can  extend  this  type  of  analysis  to 
a  real  data  set. 
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4  Results 

4.1  Improved  Target  Detection  (Task  1) 

4.1.1  Introduction  and  Background 

The  objective  of  this  work  was  to  improve  upon  the  process  of  mapping  and  picking  targets  for 
cued  re-acquisition.  The  method  was  proposed  and  developed  because  authors  believed  there  is 
a  better  way  to  pick  and  select  targets  for  cued  re-acquisition. 

The  fundamental  input  data  for  the  process  is  dynamic  survey  data  that  is  collected  along  survey 
lines.  This  data  is  typically  different  than  static  data  used  for  target  classification  in  that  it 
typically  is  for  just  one  transmitter  loop  (the  Z-axis  loop),  it  has  short  data-point  collection  times 
and  short  decay-transient  times.  The  data  collection  intervals  for  each  data  point  are  typically 
0.1s  and  the  signal  repetition  rates  are  typically  90  or  270  Hz  (2.7ms  or  0.9ms  transient  decay 
times). 

The  objective  of  analysis  of  the  data  is  to  produce  a  list  of  targets,  cued  for  re-acquisition.  The 
process  involves  selection  of  spatial  anomalies  in  the  data  and  production  of  the  cued-target  list. 
There  is  usually  no  shortage  of  spatial  anomalies  that  can  be  detected.  But  re-acquisition  of  all 
of  the  anomalies  would  be  cost-prohibitive.  Therefore,  typically  a  threshold  is  included  in  the 
process  to  select  cued  targets  from  a  map  of  all  possible  targets.  The  threshold  can  be  set  by 
ancillary  measurements  to  assure  selection  of  all  anomalies  for  the  ‘weakest’  target  of  interest, 
based  on  historical  evidence  of  site  usage  or  exploratory  surveys. 


4.1.2  Traditional  method 

Simplified,  the  traditional  method  plots  signal  amplitude  on  a  2D  map  and  picks  peaks  of  signal 
level. 

•  Sequential  data  points  are  collected  as  the  system  is  moved  across  the  ground.  A  data 
point  consists  of  one  or  more  signal  amplitudes. 

•  The  signal  amplitude  is  plotted  at  a  point  on  a  map.  For  traditional  EM61  data,  signal 
amplitude  is  plotted  at  the  center  point  of  the  loop  at  the  time  that  each  data  point  is 
collected.  Dynamic  surveys  with  MetalMapper  have  been  similar,  as  described  below. 

•  The  data  are  processed  by  edits,  filtering,  background  subtraction,  coordinate  analysis 
and  orientation  corrections,  and  other  considerations. 

•  The  data  points  are  then  plotted  on  a  2D  map. 

•  The  map  is  ‘gridded’  and  ‘contoured.’  This  step  and  the  previous  two  steps  are  often 
iterative. 

•  Peaks  are  picked  from  the  2-D  contour  plot.  Picking  of  peaks  usually  involves  selection 
of  an  amplitude  threshold  -  only  peaks  greater  than  the  threshold  are  selected.  The 
targets’  X-Y  locations  are  taken  to  be  the  coordinates  of  the  peaks  and  the  target 
amplitudes  are  taken  to  be  the  amplitudes  of  the  peaks. 
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4.1.3  MetalMapper  refinements  to  traditional  method 

MetalMappers  were  used  to  collect  dynamic  data  over  the  San  Luis  Obispo  and  Camp  Butner 
sites. 

When  ‘traditional  method’  maps  are  made  from  MetalMapper  data,  the  method  just  discussed  is 
slightly  more  complicated  because  the  MM  has  seven  small  sensors  that  sense  the  magnetic  field 
at  seven  separate  points  within  the  Z-transmitter  loop.  For  shallow  targets,  as  the  MM  array  is 
moved  along  a  line,  the  signal  amplitude  from  each  sensor  reaches  a  peak  when  the  sensor  is 
nearest  the  target5.  On  a  map,  when  the  sensor’s  signal  amplitude  is  plotted  at  the  coordinate  of 
the  center  of  the  antenna  array,  the  sharp,  narrow  peaks  of  signal  amplitude  for  each  sensor  are 
spatially  offset  with  respect  to  one  another  because  the  sensors  in  front  reach  the  target  first.  We 
would  rather  that  the  peaks  all  line  up  with  one  another  since  they  are  each  sensing  a  single 
target.  Contrarily,  for  a  deep  target,  the  peaks  from  all  sensors  are  broader  and  appear  nearly 
aligned. 

Because  of  this,  for  MM  data,  we  have  historically  used  either  of  two  methods  to  produce  a  map 
for  target  detection  and  selection. 

1.  The  5IZ  method:  We  used  the  five  innermost  sensors  and  computed  the  average  signal 
from  the  Z  coils  of  those  sensors.  This  value  is  plotted  at  the  center  of  the  antenna 
array  since  the  receivers  are  symmetrically  distributed  about  the  center.  This  method 
solves  the  peak-offset  problem  by  producing  an  average  signal  that  is  equivalently 
from  a  larger  loop  and  that  has  a  broader  peak.  This  method  approximates  the  signal 
that  could  be  received  using  an  EM61  because  we  are  using  only  the  Z  component, 
like  an  EM61  sensor,  and  because  we  are  averaging  the  signal  over  a  spatial  area  more 
nearly  the  size  of  an  EM61  receiver  loop. 

2.  The  ‘Split  Cube’  method:  In  this  method,  the  data  from  each  cube  is  plotted  at  the 
coordinate  where  the  sensor  is  located  at  the  time  of  the  data  point.  Since  each  cube 
follows  its  own  track  as  the  antenna  array  is  moved  across  the  ground,  this  method 
essentially  ‘splits’  a  single  survey  line  into  seven  survey  lines.  The  data  from  the 
seven  lines  are  contoured  on  a  2D  map  as  usual.  For  shallow  targets  this  method 
works  reasonably  well  because  it  aligns  the  peak  offsets  along  the  profile  as  discussed 
just  previously.  However,  for  deeper  targets,  a  sensor’s  peak  response  may  be  nearer 
to  the  center  of  the  antenna  array  than  to  the  coordinate  of  the  sensor.  So,  one 
refinement  to  this  method  has  been  to  use  an  effective  coordinate  for  each  sensor  that 
is  part  way  between  the  center  of  the  antenna  array  and  the  sensor  itself.  Within  this 
method  we  have  used  either  the  Z  signal  from  each  sensor  or  the  magnitude  (M)  signal 
from  each  of  the  seven  sensors,  where  magnitude  is  the  orthogonal  summation  of  the 
Z,  Y,  and  X  components 


5  This  is  not  precisely  true  because  the  target  response  is  a  function  of  the  amplitude  of  the  primary  field.  For 
shallow  targets  the  primary  field  amplitude  is  largest  near  the  edge  of  the  transmitter  loop.  Thus  signal  amplitude  at 
a  given  sensor  is  a  function  of  the  distance  of  the  target  from  the  edge  of  the  transmitter  loop  and  the  distance  of  the 
sensor  from  the  target. 
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4.1.4  Objective  of  this  work  and  basic  considerations 

In  this  project  we  attempt  to  improve  upon  the  traditional  method  of  mapping  and  target  picking. 
We  believed  from  observation  of  significant  amounts  of  data  improved  maps  could  be  made.  For 
example,  it  stands  to  reason  that  a  small  shallow  target  can  be  identified  because  it  shows  a 
signal  anomaly  in  only  two  or  three  sensors  that  pass  near  it,  while  a  deep  target  can  be  identified 
because  it  shows  a  signal  anomaly  of  roughly  equal  magnitude  in  all  seven  sensors.  In  both 
cases,  the  signal  amplitude  might  be  about  the  same;  but  we  know  that  the  deeper  target  has  a 
substantially  larger  moment,  and  that  this  moment  becomes  amplified  even  more  if  we  consider 
the  size  of  the  transmitted  field  that  activated  the  target. 

The  traditional  method  of  producing  a  cued-target  list  can  be  thought  of  as  a  two-step  process.  In 
the  first  step,  a  data  point  is  calculated  from  each  measurement  point  and  assigned  a  coordinate 
on  the  map  -  this  data  point  is  not  yet  a  pick  because  no  criteria  are  applied.  Traditionally,  this 
data  point  is  assigned  the  coordinate  of  the  center  of  the  antenna  array  and  the  amplitude  of  the 
received  signal.  The  results  from  all  of  the  data  points  are  plotted  and  the  map  is  gridded  and 
contoured.  Finally,  targets  are  picked  from  peak  amplitudes  of  anomalies  on  the  2D  map6.  This 
second  step  is  the  picking  process  where  the  criterion  is  target  amplitude. 

The  5IZ  method  for  the  MetalMapper  is  an  identical  process  where  the  amplitude  of  the  target  is 
the  taken  as  the  average  Z-component  for  the  five  innermost  sensors.  The  ‘Split  Cube’  method  is 
similar  except  that  seven  targets  are  picked  from  the  data.  Each  of  the  targets  is  picked  at  the 
location  of  an  individual  sensor  and  each  is  assigned  the  signal  amplitude  at  that  sensor.  In  either 
case,  previous  processing  of  MetalMapper  dynamic  data  has  been  not  dissimilar  to  traditional 
methods,  nor  has  it  been  a  substantial  improvement. 

4.1.5  Enabling  and  limiting  concepts 

An  important  point  is  that  the  MetalMapper  collects  vector  data  at  each  sensor.  Yet  the  methods 
used  to  analyze  the  San  Luis  Obispo  and  Camp  Butner  dynamic  surveys  did  not  make  use  of  the 
MetalMapper’s  vector  capabilities.  Following  the  reasoning  above,  our  initial  objective  is  to 
pick  zero  or  more  tentative  targets  from  each  data  point  and  to  assign  each  target  a  coordinate 
and  a  scalar  magnitude.  Our  final  objective  is  to  examine  the  map  of  scalar  magnitudes  and  to 
pick  cued  targets  for  reacquisition.  So  to  take  advantage  of  the  MetalMapper’s  vector 
capabilities,  somewhere  within  the  method  we  must  produce  one  or  more  scalars  from  one  or 
more  vectors.  The  work  herein  has  used  physics-based  approaches  to  use  the  vector  data 
produced  by  the  MetalMapper.  This  is  the  fundamental  enabling  concept  of  the  proposed 
method. 

This  work  has  been  limited  to  following  a  procedure  analogous  to  the  conventional  two-step 
process.  That  is,  we  are  considering  data  points  along  a  line  independently  of  one  another  in  the 
first  step,  and  considering  all  data  points  together  in  the  second  step.  The  proposed  method 
provides  a  substantial  improvement  to  the  first  step  and  changes  the  methodology  of  the  second 
step.  In  the  first  step  we  select  zero  or  more  targets  from  each  data  point  and  assign  each  target  a 
coordinate  and  amplitude.  But  the  selections  are  made  (using  a  minimum  correlation  factor 


6  This  of  course  is  a  grossly  simplified  description  since  editing,  subtracting  backgrounds,  smoothing  and  filtering, 
companding,  and  plotting  are  usually  often  complicated  activities.  In  addition,  the  process  of  picking  peaks  from  the 
2D  map  can  become  very  complicated. 
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criterion),  they  are  assigned  an  estimated  coordinate,  and  their  amplitude  is  a  function  of  their 
depth  and  the  magnitude  of  the  primary  energizing  field  at  that  coordinate.  It  turns  out  that  the 
picks  tend  to  group  or  cluster  so  well  that  a  gridded  contour  map  of  those  picks  is  not  helpful.  So 
instead,  the  second  step  is  one  of  applying  a  clustering  algorithm  to  produce  the  final  set  of  cued 
targets,  instead  of  gridding  and  peak-picking. 

The  number  of  tentative  targets  from  an  individual  data  point  is  an  important  factor.  A 
traditional  survey  produces  exactly  one  data  point  from  each  measurement  point.  Our  traditional 
split  cube  method  produces  seven  data  points  from  each  measurement  point.  The  method 
proposed  herein  produces  zero  or  more  targets  from  each  data  point.  The  number  of  possible 
targets  is  fundamentally  related  to  the  diversity  of  the  data.  The  MetalMapper  measures  seven  3- 
component  signals  —  21  components  that  are  not  completely  independent.  So  an  important 
question  is  How  many  targets  can  be  effectively  separated  and  detected  from  one  data  point? 

Given  uncertainties  and  noise  in  the  data,  and  given  our  experience  in  experiments  in 
performance  of  this  work,  we  believe  the  practical  answer  is  one  or  occasionally  two  or  rarely 
three  for  the  shallowest  targets,  and  only  one  for  any  deeper  target.  If  there  are  interfering 
targets,  the  proposed  method  may  be  able  to  detect  that  fact  but  cannot  separate  the  multiple 
targets,  if  for  no  other  reason  than  the  lack  of  diversity  in  the  data. 

4.1.6  Requirements  placed  upon  our  methods 

In  the  method  we  have  chosen,  we  divide  a  volume  below  the  sensor  array  into  a  set  of  voxels. 
Once  we  do  this,  our  next  task  is  to  produce  a  statistic  for  each  voxel  that  relates  to  the  likelihood 
that  the  voxel  contains  a  target.  Then  our  final  task  is  to  choose  which  voxels  are  most  likely  to 
contain  a  target  and  to  assign  values  to  those  choices  that  are  related  to  the  moment  of  the  target. 
Note  that  this  method  is  a  3D  method,  as  differentiated  from  the  conventional  2D  method.  This 
complicates  choices  for  data  visualization.  For  most  of  this  work  we  have  generated  groups  of 
plots  to  visualize  the  data.  The  plots  are  shown  as  horizontal  2D  maps  for  each  layer  of  voxels. 

We  would  like  a  method  that  produces  maps  with  the  following  characteristics: 

1.  A  plot  of  some  statistically  constant  value  when  there  are  no  targets  present.  The 
statistically  constant  value  is  related  to  the  noise  level  (assumed  statistically  constant) 
from  each  of  the  sensors  and  we  assume  that  this  noise  level  is  statistically  the  same  for 
all  of  the  sensors.  The  plotted  value  should  neither  increase  with  depth  nor  increase  with 
distance  from  the  sensors.  It  should  show  a  ‘background’  map  that  is  random  but  low 
level,  like  the  noise  level  from  the  sensors.  Since  a  map  for  a  given  data  point  is 
produced  from  the  data  for  only  seven  sensors,  an  individual  map  cannot  show 
randomness  very  well.  The  randomness  must  be  perceived  from  sequences  of  maps. 

2.  A  plot  that  shows  an  increase  of  the  plotted  value  when  signals  from  one  or  more  sensors 
rise  above  background  levels.  We  would  hope  that  the  plotted  value  shows  some 
indication  of  the  target’s  moment.  Because  a  target’s  moment  increases  by  R-cubed  for  a 
constant  observed  signal,  a  better  requirement  is  that  the  mapped  value  be  monotonically 
related,  but  not  necessarily  proportional,  to  the  observed  signal. 

3.  A  plot  that  shows  an  increase  of  the  plotted  value  for  voxels  that  are  more  likely  to 
contain  a  target.  This  likelihood  is  a  function  of  the  agreement  between  observed  signals 
from  two  or  more  sensors.  If  the  observed  signals  from  two  or  more  sensors  correlate 
well  with  the  signals  that  would  be  produced  from  a  target  in  a  cell  under  consideration, 
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then  a  larger  value  should  be  plotted.  These  considerations  must  be  true  both  laterally 
and  vertically,  because  we  assume  that  we  are  going  to  pick  one  or  more  targets  from  the 
maps. 

4.  A  smaller  anomalous  area  of  increased  values  when  the  target  is  shallow  and  a  larger 
range  of  increased  values  when  the  target  is  deep.  This  correlates  to  the  idea  that  a 
shallow  target  can  have  a  signal  that  couples  into  only  two  or  three  sensors,  whereas,  a 
signal  from  a  deep  target  will  couple  into  most  of  the  sensors. 

4.1.7  Progress  of  our  work  and  lessons  learned 

The  requirements  we  just  stated  were  partly  a  result  of  our  work  as  it  progressed.  In  our  studies 
we  pursued  the  following  methods. 

4.1.7. 1  Sensor  groups  and  forward  computations 

We  originally  proposed  to  use  pairs  and/or  groups  of  sensors  to  compute  an  apparent  target. 
That  computation  was  to  be  followed  by  examination  of  apparent  targets  to  determine  if  any 
them  were  likely  to  be  real  and  if  so,  how  big  they  were.  This  approach  was  based  on  idea  that 
fundamentally  two  sensors  are  capable  of  detecting  a  target,  its  position,  and  it’s  moment  -  six 
measurements  (Bix,  Biy,  Biz,  B2X,  B2y,  B2Z)  and  six  unknowns  (x,  y,  z,  Mx,  My,  Mz).  The  original 
concept  was  that  shallow  targets  could  be  detectable  by  just  two  sensors  whereas  deep  target 
would  require  all  sensors. 

This  method  was  problematic  in  terms  of  defining  and  using  sensor  groups.  We  experimented 
with  using  two  specified  sensors  for  specific  triangular  areas  in  shallow  layers,  and  all  sensors 
for  the  whole  lateral  area  for  deep  layers,  and  other  combinations  for  medium  depth  targets.  But 
given  that  a  target’s  depth  is  unknown  at  the  outset,  selection  of  sensor  groups  is  not  straight 
forward. 

Furthermore,  we  were  unable  to  develop  an  analytical  method  for  the  solution  to  the  two-sensor 
problem.  And  the  three-or-more  cube  problem  is  over-determined  and  requires  some  method  of 
minimization. 

Looking  at  sensor  groups  convinced  us  that  a  better  approach  was  to  use  all  sensors,  all  of  the 
time,  but  to  weight  the  contribution  of  each  sensor  in  computations  so  that  closest  sensors 
contributed  most  to  any  result.  Given  that  we  were  performing  computations  for  a  given  voxel, 
we  could  compute  the  distance  to  each  sensor  and  weight  its  contribution  to  a  result  as  a  function 
of  its  distance.  Given  that  a  magnetic  dipole  field  decreases  as  1/R3  it  seemed  intuitive  to 
establish  weighting  as  1/R3.  That  way,  for  shallow  targets,  the  closest  one  or  two  cubes  would 
dominate  the  solution  but  for  deep  targets  all  sensors  would  contribute  to  the  solution. 
Furthermore  this  approach  allows  a  kind  of  continuous  transition  between  use  of  just  a  few 
sensors  and  use  of  all  sensors. 

4.1.7. 2  Method  of  Apparent  Moments 

This  second  approach  and  the  approach  that  required  a  substantial  amount  of  effort  on  this 
project  we  called  the  method  of  Apparent  Moments.  However,  this  method  proved  to  be  severely 
limited  and  was  not  implemented  in  the  final  analysis.  It  is  described  briefly  here. 

Our  overall  objective  in  any  method  was  fundamentally  to  determine  position  and  apparent 
moment  of  a  target  (or  targets)  given  one  measurement  of  MetalMapper  data.  Assuming  a  dipole 
target,  and  assuming  a  known  position  of  that  dipole,  computations  of  its  magnetic  field  are 
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simple  and  straight  forward.  Importantly,  given  the  position  of  target,  and  given  its  magnetic 
field,  its  moment  (magnitude  and  orientation)  is  a  simple  forward  computation. 

Given  this,  with  the  volume  below  the  array  divided  into  a  grid  of  voxels,  our  approach  was  to 
compute  an  apparent  moment  in  each  voxel  for  each  signal  observed  at  the  seven  sensors.  Given 
seven  sensors,  the  result  was  seven  apparent  moments  in  each  voxel.  From  here,  the  concept  was 
to  examine  the  apparent  moments  in  each  voxel.  A  good  target  would  produce  seven  moments 
all  having  with  same  orientation  and  magnitude  whereas  noise  would  produce  seven  moments 
with  random  orientations  and  magnitudes.  So  the  objective  became  one  of  developing  a  method 
to  determine  whether  the  seven  apparent  moments  matched  each  other  or  whether  they  did  not. 

We  computed  a  mean  apparent  moment  in  each  voxel  as  a  weighted  average  of  the  seven 
apparent  moments.  Weighting  was  chosen  as  discussed  above  to  amplify  the  contribution  from 
closest  sensors  to  a  voxel  under  consideration.  This  produces  an  apparent  mean  moment  in  each 
voxel  but  it  does  not  say  which  mean  moment  is  highest  quality.  If  the  largest  moment  is 
selected,  the  deepest  voxels  are  almost  always  chosen  because  they  have  the  largest  mean 
moments  even  though  the  individual  moments  do  not  agree.  To  select  a  ‘best’  mean  moment 
requires  computation  of  a  quality  factor. 

We  developed  a  quality  factor  based  on  the  correlation/matching  of  individual  apparent 
individual  moments  in  a  voxel.  We  optionally  included  another  quality  factor  based  upon 
whether  or  not  the  apparent  moment(s)  in  a  voxel  had  an  expected  orientation  -  i.e.  whether  or 
not  they  matched  the  orientation  of  the  primary  field  for  that  voxel.  We  experimented  with 
methods  to  compute  these  quality  factors  and  to  select  the  voxel(s)  with  the  best  quality. 

None  of  the  methods  we  tried  performed  particularly  well.  A  constant  difficulty  with  the  method 
was  that  dipole  moments  computed  for  deep  voxels  were  (not  unexpectedly)  always  larger  than 
moments  computed  for  shallow  voxels.  It  was  difficult  to  compute  quality  factors  that  robustly 
picked  targets.  We  found  that  the  signal-pattern  correlation  method  in  the  next  section 
consistently  outperformed  the  apparent  moment  method.  Nevertheless,  the  apparent  moment 
method  allows  a  user  to  pick  targets  without  assuming  that  the  target  dipole  is  aligned  with  the 
primary  field.  If  this  issue  becomes  more  important  than  believed,  the  apparent-moment  method 
should  be  re-visited. 

We  obtained  the  best  results  with  a  method  that  suggested  itself  during  this  work  and  is  described 
next. 

4.1.8  Method  of  Signal  Pattern  Correlation 

We  call  this  method  Signal  Pattern  Correlation  because  it  employs  a  simple  correlation  of  the 
observed  signal  vectors  with  a  set  of  reference  signal  vectors.  The  method  might  be  called 
Pattern  Matching  or  even  Matched  Filtering. 

4.1.8.1  Computations 

This  method  is  based  upon  a  correlation  of  a  set  of  reference  signals  with  the  observed  signals 
for  each  of  the  seven  sensors.  We  divide  the  volume  below  the  antenna  array  into  a  set  of  cells 
or  voxels  and  then  perform  a  correlation  for  every  voxel  in  our  test  volume.  We  then  choose 
zero  or  more  voxels  having  highest  correlations  and/or  largest  signal  levels. 

We  assume  the  traditional  dipole  model  and  compute  the  signal  that  would  be  observed  at  each 
sensor  from 
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Sy kl  ~  n3  fsijkl^ijk  ^  ~  1’^,. 


ijk 


where  Sijkl  is  the  signal  at  sensor  /  from  an  object  at  location  (i,j,k),  and 

2  —  3r  rT  —  I 

oijkl  -}lijkllijkl  1 


(4.1) 


(4.2) 

where  ijk  are  the  indices  of  cell  containing  the  target,  R  is  the  vector  from  the  cell  at  ijk  to  the 
sensor  l,  R  is  the  magnitude  of  R,  r  is  the  unit  vector  aligned  with  R,  M///{  is  an  arbitrary  moment 

representing  the  target,  and  M(.//{  is  shown  with  indices  to  indicate  that  it  is  located  in  cell  ijk. 
We  assume  that  a  target  of  interest  will  have  the  same  orientation  as  the  primary  field.  Thus 

M*  =  Mlijk  (4.3) 

where  Tijk  is  the  primary  transmitter  field  at  the  cell  ijk,  T  is  its  magnitude,  and  M  is  a  constant 

representing  the  magnitude  of  the  target.  The  assumption  that  the  target  moment  has  the  same 
orientation  as  the  primary  field  is  correct  in  early  time  channels  for  any  target  and  correct  for  all 
time  channels  for  a  spherical  target. 

To  compute  reference  vectors  to  be  used  for  correlation,  we  set  M  =  1  and  compute  a  reference 
signal  Sijkl  from  Equation  4.1.  Then  we  correlate  the  observed  signal  S;  with  Sijkl  We  chose  a 

dot  product  as  the  method  of  correlation.  We  note  that  a  dot  product  of  seven  3 -component 
vectors  is  the  same  as  a  simple  correlation  of  21  signals  with  21  references.  The  correlation 
factor  for  a  given  voxel  is 
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where  Sm  is  the  reference  vector  for  cell  ijk  and  sensor  /, 
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Note  that  a  s  is  like  the  root-mean-square  (RMS)  value7  of  the  signals  at  seven  sensors  and  that 
ojjk  is  a  normalizing  factor  to  cause  the  maximum  correlation  to  be  unity. 


If  an  actual  target  is  in  the  cell  ijk,  then 

S;  =  MSyu  5  Cy kl  =  ^  5  and  M  =  ■ 

This  equation  gives  the  moment  of  a  dipole  if  it  is  actually  located  in  cell  ijk. 


(4.8) 


7  Actually,  is  a  ‘root-sum-square’  (RSS)  value.  In  this  case  case  RMS  and  RSS  are  related  by  the  constant  factor 
v  7  for  seven  sensors. 
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Note  that  because  the  reference  signals  are  scaled  by  1  /R3  in  Equation  (4.1),  the  weighting 
discussed  earlier,  to  use  only  nearest  sensors  for  shallow  targets,  and  to  use  all  sensors  for  deep 
targets,  is  implicitly  included. 

4.1.8.2  Selection  of  Targets 

Once  the  values  above  are  computed  for  each  voxel,  the  next  step  is  to  select  the  ‘best’  voxel(s) 
as  the  one(s)  most  likely  to  contain  a  target.  For  this  we  examine  the  correlation  coefficient 
knowing  that  a  perfect  correlation  is  +1.  Also  note  that  for  a  single  data  point,  we  have  a  single 
realization  of  the  received  signals,  and  that  the  effective  amplitude  of  the  received  signals  is  as. 
Thus  we  have  two  parameters  that  can  be  used  to  decide  whether  target(s)  are  present.  The  first, 
as  is  an  enabling  indicator  -  if  the  signal  amplitude  is  small,  such  as  ‘background’  then  we  know 
there  are  no  targets  present  even  if  noise  should  indicate  a  high  correlation  coefficient  in  any 
voxel.  If  the  signal  amplitude  is  large  enough,  such  as  above  a  threshold,  then  we  know  to 
examine  the  correlation  coefficients  in  all  of  the  voxels  to  determine  where  the  target(s)  might 
be. 

One  method  is  to  select  the  best  voxel  in  the  entire  volume  of  voxels  -  the  voxel  that  had  the 
largest  correlation  coefficient.  This  method  is  modified  by  selecting  no  voxels  if  the  best  voxel 
was  on  the  edge  (or  surface)  of  the  voxel  array  -  the  justification  is  that  if  the  correlation  was 
spatially  increasing  toward  any  edge  or  surface  of  the  array,  then  the  maximum  voxel  might  be 
located  outside  the  array  -  therefore  the  target  is  not  located  within  the  voxel  array.  For 
discussion  purposes  herein  we  have  named  this  method  the  best-pick  method. 

A  second  method  is  to  select  any  voxel  having  a  spatially  local  maximum  correlation  coefficient. 
Any  voxel  under  consideration  has  26  neighbors  (a  3x3x3  array  in  which  the  voxel  of  interest  is 
in  the  center).  Therefore  the  method  requires  that  all  26  neighbors  have  a  correlation  coefficient 
less  than  the  center  voxel.  A  characteristic  of  this  method  is  that  it  is  capable  of  selecting 
multiple  targets  from  the  data.  We  have  named  this  method  the  multi-pick  method. 

Both  of  the  methods  can  be  modified  by  requiring  not  only  that  the  signal  level  £TV  be  above  a 
threshold,  but  also  that  the  correlation  coefficient  c,  .*  be  above  a  threshold.  These  thresholds  are 
a  key  characteristic  of  the  method.  They  allow  completely  automated  selection  of  targets  with 
deterministic  specifications  selection  parameters,  similar  to  the  conventional  specification  of  a 
threshold  when  picking  peaks  from  a  contour  map. 

Once  a  target  is  selected,  its  moment  is  computed  from  Equation  4.8.  Note  that  selection  of 
targets  depends  only  on  signal  amplitude  and  on  correlation  of  the  signal  pattern  to  an  expected 
pattern.  In  this  regard  it  is  not  completely  dissimilar  to  the  traditional  method  of  signal  contour 
mapping  and  picking  peaks.  However,  this  method  adds  and  allows  an  important  refinement  to 
the  traditional  process.  Once  a  target  has  been  picked,  its  moment  is  determined.  So  a  second 
level  of  target  picking  can  be  done  by  examining  moments  of  the  targets.  The  conventional 
method  never  produces  the  moment  of  the  targets  nor  does  it  estimate  target  depth.  This  new 
method  allows  the  analyst  to  at  first  pick  all  targets  having  adequate  signal  amplitude  at  the 
sensors,  and  then  to  substantially  reduce  the  number  of  picks  based  on  target  moment,  instead  of 
signal  amplitude.  The  second  step  improves  on  the  issue  of  picking  too  many  near-surface 
targets  at  the  expense  of  missing  deeper  targets. 

This  method  tends  to  display  points  on  a  map  that  cluster  into  groups  because  a  given  target  can 
be  detected  in  many  individual  data  points.  Whereas  in  the  traditional  method  all  of  the 
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detections  are  smeared  along  a  survey  line,  in  this  method  the  detections  tend  to  cluster  in  groups 
directly  at  the  target.  We  found  that  if  we  allowed  ‘detections’  based  on  small,  noise-only 
signals,  the  results  tended  to  scatter  as  expected.  But  when  we  plotted  the  resulting  moments  on 
a  map,  the  (scattered)  deep  targets  dominated  the  map  because  of  their  larger  amplitudes. 
Therefore  we  employed  a  signal  threshold  discussed  earlier.  This  causes  many  areas  on  the  map 
to  have  no  data  points.  The  maps  tend  to  have  spatially  small  clusters  of  points  directly  where  a 
target  is  located  and  not  many  points  elsewhere.  Thus  the  traditional  method  of  gridding  and 
contouring  is  not  very  helpful.  Instead,  a  method  is  needed  to  pick  clusters  of  points  and  to 
determine  the  mean  location  and  value  (; moment  in  this  case)  of  each  cluster.  Of  course  a 
clustering  algorithm  becomes  complex  to  determine  where  there  are  clusters  and  which  points 
belong  to  each  cluster,  in  the  same  way  that  picking  peaks  from  a  2D  map  becomes  complex  to 
determine  which  maxima  are  independent  peaks.  A  simple  clustering  algorithm  was 
implemented  for  this  study  but  is  not  part  of  the  proposed  method.  If  future  work  is  done  with 
this  method,  selection  and  implementation  of  better  clustering  algorithms  will  be  part  of  that 
work. 

4.1.8.3  Implementation  and  Performance 

The  method  was  implemented  in  two  ways.  The  second  method  below  is  the  most  important 
because  its  result  is  to  improve  the  process  of  producing  a  cued  target  list.  However,  the  method 
is  also  useful  for  user/operator  feedback,  in  real  time,  in  the  field.  That  implementation  is 
presented  first. 

4.1.8.4  Real-time  detection  and  dancing  arrows  display  in  program  EM3D. 

The  MetalMapper  data  acquisition  program  EM3D  has  traditionally  provide  a  display  known  as 
dancing  arrows.  This  display  has  been  used  to  center  the  MetalMapper  over  a  target  because  the 
arrows  tend  to  point  directly  toward  the  target.  The  operator  can  also  see  targets  passing  by  the 
array  through  observation  of  the  dancing  arrows. 

The  correlation  detection  method  in  this  report  was  implemented  in  EM3D  to  add  detected 
targets  to  the  dancing  arrows  display.  The  method  selects  zero  or  more  criteria  discussed  above 
and  plots  the  target(s)  along  with  the  dancing  arrows.  The  detected  depth  is  displayed  as  a  color 
and  the  magnitude  of  the  moment  of  the  target  is  displayed  as  symbol  size.  The  result  is  that  the 
user  can  watch  targets  pass  by  the  array,  and  will  have  estimates  of  their  size  and  depth. 

The  method  is  implemented  as  an  array  of  voxels  that  is  1.56m  wide  with  cells  that  are  6.5cm 
square  and  20cm  tall.  The  center  of  the  first  cell  is  positioned  5  cm  above  the  center  of  the  Z- 
antenna  and  receiver  array  so  that  the  first  layer  of  cells  not  on  the  surface  of  the  voxel  volume  is 
15cm  below  the  center  of  the  Z  antenna,  or  roughly  ground  level.  The  voxel  array  contains  25  x 
25  x  7  cells.  The  effective  volume  of  useful  cells  is  1.43m  x  1.43m  x  100cm,  from  15cm  to  1 15 
cm  below  the  antenna  array.  The  entire  array  contains  4375  cells  and  can  be  run  for  every  data 
point  collected  by  EM3D,  without  slowing  down  data  collection,  when  data  points  are  collected 
at  0.1s  intervals.  Tests  were  not  made  to  determine  how  much  faster  than  this  the  routine  runs. 

Samples  of  screen  shots  from  EM3D  are  shown  in  Figure  4.1.  It  shows  the  dancing  arrows 
panel.  The  pictures  are  shown  from  real  data  collected  along  Line  1,  the  westernmost  line,  of 
SLO.  The  target  in  the  left  hand  panel  is  one  of  the  smaller  targets  near  the  south  end  of  the  line. 
The  target  in  the  right  hand  panel  is  the  larger  anomaly  about  mid-way  and  to  the  east  of  the  line. 


22 


Expanded  Processing  Techniques  for  EMI  Systems 


SERDP  MR- 1772 


Figure  4.1:  Screen  shot  pictures  of  implementation  in  EM3D.  The  left  panel  shows  a  target  at  a  depth  of  35cm  passing 
under  the  array.  The  right  hand  panel  shows  a  much  larger  target  at  a  depth  of  75cm  or  so  passing  just  outside  the  array. 


4.1.8. 5  Mapping  detection  and  display  in  Oasis  montaj 

The  method  has  been  implemented  as  a  GX-Net  add-on  routine  for  Geosoft’s  Oasis  montaj.  The 
implementation  contains  both  the  best-pick  and  multi-pick  methods. 

In  Oasis  montaj,  it  is  convenient  to  show  the  best  pick  for  each  data  point  and  to  display  that  data 
beside  the  rest  of  the  data  for  each  fiducial.  Displaying  multiple  picks  per  data-fiducial  requires 
multiple  channels  or  an  array  channel  with  variable  length.  Therefore  multiple  picks  are  entered 
in  a  new  line  or  new  group  in  Oasis  montaj.  The  user  must  select  either  original  survey  lines  or 
the  one  multi-pick  line/group  when  plotting  picks  on  a  map. 

The  method  allows  the  user  to  specify  the  overall  size  and  cell  size  of  the  voxel  array;  and  it 
allows  the  user  to  specify  thresholds.  The  user  interface  is  shown  in  Figure  4.2. 

The  implementation  returns  correlation  coefficient,  coordinate  offsets  and  moment  for  the 
target(s)  selected.  Within  montaj,  the  coordinate  offsets,  which  are  with  respect  the  cart’s 
coordinate  system,  are  manipulated  according  to  cart  heading  and  orientation,  to  produce  a 
ground-referenced  coordinate  for  each  target. 

Within  montaj,  the  user  must  compute  absolute  coordinates  for  the  target,  given  cart-relative 
coordinates  of  the  target(s)  and  orientation  of  the  cart.  For  this  study,  a  simple  correction  based 
on  cart  heading  was  used,  where  cart  heading  was  the  computed  true  heading,  using  GPS 
coordinates. 
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Figure  4.2:  User  interface  for  Signal  Pattern  Correlation  within  Oasis  montaj. 


Performance  of  the  method  is  shown  in  Figure  4.3.  Results  showing  the  intended  use  of  the 
method  are  shown  in  the  next  section  on  Results. 

Preparation  of  data  within  montaj  is  a  non-trivial  process.  Coordinates  must  be  imported  and 
converted  from  Longitude/Latitude  to  UTM,  data  must  be  edited,  orientation  data  must  be 
filtered,  and  EM  data  must  be  gated  or  other-wise  manipulated  to  produce  scalar  values  at 
individual  coordinates,  so  that  it  can  plotted  on  a  map. 

Figure  4.3  shows  application  of  the  method  along  one  line.  This  line  is  the  westernmost  line  in 
the  figures  shown  later.  The  top  panel  shows  the  traditional  5IZ  signal  in  green  and  the  root- 
sum-square  (RSS)  signal.  Although  labeled  ‘SignalRMS’  and  ‘Avg5IZ,’  these  signals  are 
actually  the  sum  of  5  inner  Z  signals  and  the  root-sum-squares  of  the  seven  sensor  signals.8 


8  Work  on  this  method  leads  the  author  to  wonder  why  we  have  never  (traditionally)  used  the  RMS  or  RSS  signal  to 
produce  the  contour  maps  for  conventional  peak  picking.  It  appears  to  be  a  reasonably  robust  signal  that  should 
produce  a  better  map  than  the  methods  we  have  previously  used. 
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Thus,  the  orange  trace  is  the  signal,  as ,  in  Equation  (4.7).  The  other  four  panels  show 
computations  from  the  method  for  the  best  pick.  The  second  panel  shows  the  correlation 
coefficient  and  the  dipole  moment  for  the  best  pick.  For  these  plots,  we  used  an  RSS  signal 
threshold  of  1 .0  and  a  correlation  coefficient  threshold  of  0.90.  The  bottom  three  panels  show 
the  delta  coordinates,  i.e.  the  location  of  the  detected  target  coordinates  in  the  cart  reference 
system.  For  a  perfect  detection,  we  expect  to  see  a  constant  X  coordinate  that  is  the  lateral  offset 
of  the  target  on  the  cross-track  axis,  a  linearly  varying  Y  coordinate  as  the  cart  passes  the  target 
at  constant  speed,  and  a  constant  depth.  The  results  show  that  the  method  is  indeed  doing  a  good 
job  of  detecting  and  locating  the  targets. 

Investigation  of  the  correlation  picking  method  surprised  us  in  terms  of  its  lateral  effectiveness. 
It  appears  consistently  able  to  pick  targets  up  to  0.25m  outside  of  the  survey  loop,  meaning  that 
the  loop  is  covering  a  1 ,5m  swath  along  each  line.  It  will  detect  large  targets  as  much  as  0.5m 
outside  the  survey  loop.  Thus,  while  we  initially  expected  to  implement  a  voxel  array  the  same 
size  as  the  MetalMapper  array,  lm  x  lm,  we  ended  up  with  a  voxel  array  four  times  as  big,  2.0m 
x  2.0m.  We  found  it  important  to  select  a  voxel  array  that  covered  a  2m  x  2m  area  because  of 
the  MetalMappers  ability  to  locate  targets  that  pass  outside  its  nominal  footprint. 
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Figure  4.3:  Data  profiles  along  Line  1,  the  westernmost  line.  The  Avg5IZ  trace  is  the  signal  that  is  similar  to  that 
produced  by  an  EM61.  The  orange  trace  is  effective  RMS  signal  from  seven  MetalMapper  sensors.  Note  these  traces  are 
plotted  on  a  log/linear  scale.  The  red  points  in  the  second  panel  are  correlation  coefficient  while  the  blue  points  are  dipole 
moments.  The  DelXDisp,  DelYDisp,  and  DepthDisp  traces  show  the  target  positions  in  cart-relative  coordinates  (cm  from 
the  center  of  the  cart). 
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4.1.9  Results 

To  demonstrate  the  method  in  montaj,  we  chose  ten  survey  lines  from  San  Luis  Obispo.  These 
are  lines  1  through  10  of  the  area  we  designated  NWA  in  our  survey  there. 

A  contour  map  of  all  ten  lines  in  the  study  area  is  shown  in  Figure  4.4.  The  map  is  described  in 
the  figure  caption.  Importantly,  the  contour  plot  is  made  of  the  traditional  5IZ  scalar  -  the 
average  value  of  the  Z  component  of  the  five  innermost  sensors. 

Figure  4.4  shows  correlation  picks  in  white.  The  picks  clearly  show  picks  of  every  significant 
peak  that  could  be  picked  from  the  contour  map.  The  points  impressively  fall  on  or  near  the 
peaks  of  the  contour  map  but  the  clusters  are  significantly  smaller  and  indicate  that  the  re¬ 
acquisition  coordinates  would  be  better  estimates  of  target  position  than  the  traditional  re¬ 
acquisition  coordinates.  In  this  figures,  most  of  the  target  picks  are  of  substantially  smaller 
targets  than  the  threshold  that  was  actually  used  at  San  Luis  Obispo.  In  production 
implementation,  a  contour  map  would  likely  not  be  produced.  The  white  points  in  Figure  4.4. 
would  be  plotted  on  a  blank  map  for  visual  confirmation.  But  the  next  step  is  the  one  that  is  a 
clear  advantage  of  this  method.  The  clusters  of  points  are  grouped  into  single  targets.  This 
clustering  allows  the  important  step  to  use  yet  additional  thresholds  to  select  points  for  re¬ 
acquisition.  One  threshold  is  a  distance  that  is  used  to  decide  which  points  should  belong  in  a 
given  cluster.  This  threshold  could  be  related,  if  not  identical,  to  the  distance  between  targets 
that  has  traditionally  been  used  to  decide  whether  there  is  one  or  two  targets  at  a  given  location. 
A  second  and  very  important  threshold  is  a  moment  threshold.  In  the  process  of  clustering 
points,  an  average  moment  can  be  computed  for  each  cluster  of  points.  Then  this  average 
moment  can  be  compared  to  a  threshold  to  determine  if  the  cluster  should  become  a  re¬ 
acquisition  point. 

This  final  step  is  demonstrated  in  Figure  4.5.  It  shows  clusters  of  targets  that  were  picked  using 
a  distance  threshold  of  lm  and  an  arbitrary  moment  threshold.  The  moment  threshold  was 
chosen  to  make  the  map  approximately  similar  to  the  actual  points  that  were  picked  for  re¬ 
acquisition  in  this  survey  area.  The  map  shows  the  finally  picked  clusters  as  white  circles.  It 
also  shows,  as  black  triangles,  the  targets  that  were  actually  picked  for  re-acquisition  at  San  Luis 
Obispo.  Study  of  this  map  shows  that  this  new  method  picked  a  few  targets  that  were  not  picked 
using  the  traditional  method  and  it  shows  a  few  points  that  were  picked  with  the  traditional 
method  but  were  not  picked  with  this  new  method.  While  we  cannot  show  this  on  these  2D 
plots,  we  also  know  the  depth  and  size  of  the  targets  picked  using  this  method. 

We  believe,  from  working  with  these  data  sets,  that  the  picks  using  this  new  method  are  better. 
Further  work  needs  to  be  done  to  support  this  conclusion  quantitatively.  We  suggest  a  further 
analysis  could  be  done  to  compare  these  picks  and  the  original  picks  to  the  static  data  and  even  to 
the  ground-truth  data  from  SLO. 
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Figure  4.4:  Color  contour  plot  of  the  area  studied.  The  plot  is  made  from  the  “5IZ”  value  -  a  signal  roughly  similar  to  the 
response  of  an  EM61.  The  faint  dots  show  the  track  of  data  points  along  each  of  10  lines.  The  white  crosses  are  computed 
target  positions.  The  whole  area  is  shown  to  the  left.  Two  overlapping  halves  are  enlarged  to  the  right  (top  on  the  left, 
bottom  on  the  right). 
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Figure  4.5:  Map  showing  picks  using  the  Signal  Correlation  method  (white  circles)  compared  to  picks  using  the 
traditional  method  (black  triangles).  Area  shown  is  the  same  as  in  Figure  4.4. 
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4.1.10  Conclusions  and  Recommendations 

This  new  method  for  cued  target  mapping  and  selection  has  some  salient  characteristics: 

•  When  this  new  method  is  used,  MetalMapper  data  is  capable  of  detecting  and  mapping 
targets  outside  of  the  array.  The  small  amount  of  work  done  in  this  study  (only  10  survey 
lines  at  San  Luis  Obispo)  suggests  that  it  may  be  able  to  reliably  detect  and  map  targets 
over  a  path  width  of  1.5m  whereas  a  ‘normal’  path  width  is  usually  assumed  to  be  lm, 
the  width  of  the  MetalMapper  transmitting  loop.  MetalMapper  survey  lines  have  been 
traditionally  been  spaced  0.75m  on  the  basis  of  a  1.0m  path  width.  The  wider  path  width 
effective  in  this  method  could  allow  lm  line  spacing  or  perhaps  line  spacing  of  1.25m  at 
some  sites. 

•  This  method  significantly  improves  a  common  difficulty  of  deciding  whether  multiple 
sightings  of  targets  from  two  or  more  lines  are  two  targets  or  one.  Where  a  single  target 
between  two  survey  lines  might  be  conventionally  mapped  as  two  targets,  one  on  each 
line,  the  new  method  will  more  likely  map  the  two  sightings  of  the  target  as  one  target. 
Similarly,  two  targets  that  are  offset  from  two  lines  but  not  between  the  two  lines  might 
conventionally  be  mapped  as  one  target  directly  between  the  two  lines,  whereas  this  new 
method  would  correctly  separate  the  two  targets. 

•  This  initial  work  indicates  qualitatively  that  the  correlation  method  is  robust.  Since  the 
correlation  is  a  linear  process  in  terms  of  the  correlation  coefficient,  we  expect  it  to 
behave  well  in  the  presence  of  noise.  In  this  regard  it  should  be  more  robust  than  non¬ 
linear  methods. 

•  The  method  was  implemented  (in  prototype  form)  in  real-time  in  the  data  acquisition 
program  EM3D  for  MetalMapper.  It  shows  targets  moving  across  the  dancing  arrows 
display.  It  shows  target  size  (i.e.  moment,  normalized  by  primary  field  at  the  target)  via 
size  of  the  target  ‘dot’  and  it  shows  target  depth  as  color  of  the  ‘dot.’  This 
implementation  will  allow  better  interpretation  of  the  ‘dancing  arrows’  display  and 
should  improve  positioning  during  cued  surveys,  perhaps  eliminating  some  portion  of 
repeat  cued  acquisitions. 

•  Although  not  explicitly  proven  because  of  the  small  size  of  the  data  set,  we  believe  the 
method  provides  significantly  improved  cued-target  locations.  Improved  cued  locations 
should  also  reduce  number  of  repeat  cued  acquisitions. 

•  The  multiple-pick  method  can  show  two  or  more  targets  when  those  two  targets  are 
separated  enough  in  an  EM  sense.  During  simulations  in  this  work,  we  observed  that  the 
best-pick  method  tended  to  follow  one  target  for  a  while  and  then  to  follow  the  other 
target,  as  the  array  passed  over  the  targets.  In  both  cases,  this  new  method  tended  to 
indicate  the  presence  of  multiple  targets. 

•  The  method  provides  a  repeatable  process  that  we  believe  is  more  deterministic  and  less 
‘artistic’  than  the  traditional  contouring/peak-picking  method.  Gridded  contour  maps  can 
be  significantly  altered  by  different  selection  of  data  reduction  and  contouring/mapping 
parameters.  Similarly,  two-dimensional  peak  picking  from  a  gridded  contour  map  can  be 
significantly  altered  by  choice  of  a  peak  picking  algorithm  and  its  parameters.  This  new 
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method  is  not  completely  immune  from  similar  issues  but  we  believe,  based  on  the  small 
number  of  results  from  this  study,  that  computing  detection  points  from  the  processed 
data  is  less  artistic  than  choosing  parameters  for  gridding  and  contouring,  and  that 
selecting  clusters  of  detection  points,  based  on  both  point  spacing  and  moment- 
magnitudes,  is  similar  in  artistry  to  peak  detection  from  a  2d  map. 

•  Most  importantly,  the  method  allows  a  final  picking  of  targets  based  on  three  distinct 
thresholds  or  parameters,  as  differentiated  from  the  conventional  method  that  sets  one 
threshold  that  is  applied  to  a  contour  map.  The  first  is  signal  level  and/or  signal-to-noise 
ratio,  the  second  is  correlation  coefficient  (indicative  of  signal  quality),  and  the  third  is 
magnetic  moment  (indicative  of  volume  or  mass).  We  believe  this  will  allow  better 
control  of  threshold  selections  to  assure  detection  of  minimum  size  targets  of  interest 
while  not  selecting  targets  of  non-interest,  especially  small  near-surface  scrap.  It  remains 
to  be  seen  how  many  fewer  targets  this  method  will  pick  compared  to  the  conventional 
method.  A  full  analysis  of  the  Camp  Butner  data  and/or  the  San  Luis  Obispo  data  is 
indicated. 

•  The  MetalMapper  data  in  this  method  is  acquired  using  only  the  Z  transmitter  loop.  This 
means  that  a  MetalMapper  X  and  Y  transmitter  loops  could  be  removed  from  the  antenna 
array  and  the  dynamic  survey  could  be  done  using  just  the  Z  loop.  With  the  X  and  Y 
loops  removed,  the  MetalMapper  would  be  reasonably  man-deployable.. 

We  make  the  following  recommendations: 

•  Provide  additional  funding  to  allow  a  full  analysis  of  this  method  using  the  entire  sites  of 
San  Luis  Obispo  and  Camp  Butner. 

•  Assuming  results  of  additional  studies  are  encouraging,  use  the  MetalMapper  for  new 
mapping  surveys.  We  believe  that  with  suitable  modification  in  deployment 
configuration,  that  costs  can  be  nearly  as  low  as  deployment  of  an  EM-61.  And  we 
believe  that  cost-savings  in  the  rest  of  the  survey  will  more  than  offset  the  additional 
initial  costs. 

•  Implement  this  method  in  a  multi-step  process,  each  customized  for  the  site  being 
surveyed.  The  signal  thresholds  can  be  set  according  to  the  noise  levels  and/or  other 
characteristics  of  the  site. 
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Figure  4.6:  Transient  data  for  cube  3  (center  cube)  for  a  profile  of  dynamic  MM  data  acquired  over  the  IVS  at 
Butner  (July  2010). 


4.1.11  Traditional  Inversion  Using  Dynamic  Data 

In  addition,  MMTargets  supports  a  4  parameter  (x,  y,  z,  and  radius)  inversion  that  can  be  applied 
to  dynamic  MetalMapper  data  to  support  the  target  detection  phase  of  remediation  efforts.  This 
procedure  executes  quickly  enough  to  support  pseudo-real  time  processing.  The  mean  squared 
data  misfit  is  useful  for  indicating  when  a  potential  target  is  present.  Combining  this  with  the 
anomaly  amplitudes  provides  a  high  quality  indicator  that  incorporates  physics  and  reduces  the 
uncertain  threshold  issues  that  arise  with  simple  amplitude  threshold-based  target  picking.  We 
have  applied  this  mode  of  processing  to  an  example  profile  acquired  over  the  IVS  at  Camp 
Butner  (Figure  4.6).  The  method  is  functionally  equivalent  to  the  method  of  Signal  Pattern 
Correlation  described  above  except  the  technique  does  not  rely  on  tabulated  Green’s  functions. 
Our  purpose  here  is  simply  to  report  that  MMTargets  has  this  capability.  The  results  of  the 
detector  output  are  shown  in  Figure  4.7.  At  each  point  along  the  detector  profile,  MMTargets 
produces  an  estimate  of  the  relative  target  position  and  its  size.  A  simple  threshold  applied  to  the 
“Detector  Indicator”  curve  identifies  source  point  clusters  corresponding  to  the  target. 
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Detection  Indicator 


Figure  4.7:  Results  of  operating  MMTargets  on  a  profile  of  dynamic  MM  data  acquired  over  the  IVS  at  Butner.  The  upper 
profile  shows  the  “detection  indicator”  (related  to  1/MSE).  The  lower  profile  plots  the  RMS  value  of  the  data  (similar  to  the 
magenta  profile  in  Figure  4.6). 
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4.2  Optimizing  The  MetalMapper  Receiver  Array  (Task  2) 


4.2.1  Background 


4.2.1.1  The  AOL  Study 

Early  work  on  optimum  arrays  was  reported  by  Grimm  and  Sprott  [2]  and  was  supported  by 
NAVEODTECHDIV.  Grimm  extended  that  work  (unpublished)  to  consider  a  configuration  with 
3  orthogonal  transmitters.  The  important  conclusions  arising  from  Grimm  and  Sprott’ s  study 
were: 


1.  In  all  cases,  arrays  with  triaxial  (3-component)  sensors  performed  better  than  similar 
arrays  with  1 -component  sensors  where  the  performance  metric  is  a  measure  of  the  error 
in  the  estimated  target  polarizability. 


2. 
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Figure  4.8:  Three  “optimal”  arrays  evaluated  by  Grimm  for  the  AOL  system  design.  The 
annotation  “3-C”  indicates  3-component  receivers.  Arrays  A  and  B  have  a  single  horizontal 
transmitter.  Array  C  has  3  orthogonal  transmitters  and  is  therefore  suitable  for  cued  ID 
measurements  from  a  single  point. 


For  arrays  with  a  single  1  m  x  1  m  horizontal  transmitter,  a  3  x  3  grid  of  3-component 
receivers  centered  inside  the  transmitter  loop  gave  the  best  performance  (Figure  4.8A). 
An  8  receiver  array  consisting  of  the  3  x  3  grid  with  the  center  receiver  omitted 
performed  nearly  as  well  as  the  array  with  the  3  x  3  receiver  grid  (Figure  4.8B). 


3.  An  array  of  3  orthogonal  1  m  x  1  m  transmitters  with  3 -component  receivers  centered  in 
each  of  the  4  quadrants  of  the  horizontal  loop  (Figure  4.8C)  also  performed  “very  well.” 


4.  Grimm’s  study  formed  the  basis  for  specifying  the  receiver  arrays  for  the  original 
Advanced  Ordnance  Locator  (AOL),  the  follow-on  AOL2  project,  and,  subsequently,  the 
MetalMapper  project.9  Grimm’s  work  was  based  on  model  parameter  estimations  using  a 
grid  of  data  points  centered  over  the  target.  He  did  not  perform  an  extensive  study  on  the 
performance  of  the  more  complex  3-transmitter  array  (Figure  4.8C)  using  data  from  a 
single  site.  Grimm  did  not  consider  the  deleterious  effect  of  uncertainties  in  relative 
position  and  attitude  errors  within  the  measurement  grid.  The  effect  of  these  errors  and 
the  need  for  extremely  precise  locations  (0[lcm])  has  been  noted  for  some  years  [3], 


9  The  AOL  project  was  awarded  to  Blackhawk  Geoservices  in  2003  by  NAVEODTECHDIV  (Contract  #N00174- 
03-C-0006).  The  follow-on  contract  was  awarded  in  2005  to  G&G  Sciences,  Inc.  The  MetalMapper  project  was 
funded  by  ESTCP  and  was  awarded  to  Geometries  in  2006. 


34 


Expanded  Processing  Techniques  for  EMI  Systems 


SERDP  MR- 1772 


4.2. 1.2  The  BUD  Study 

In  a  study  leading  up  to  the  configuration  of  the  BUD  system,  Smith  et  al  [4]  conducted  a  study 
similar  to  that  of  Grimm.  The  study  focused  on  optimizing  the  positions  of  4,  5,  and  6  single¬ 
component  receivers  using  a  sphere  model.  Their  study  considered  arrays  with  1,  2,  and  3 
transmitter  loops.  The  study  found  that  the  sensitivity  of  changes  in  receiver  position  is  low 
when  they  used  a  spherical  target.  However,  when  highly  elongated  (1000: 1  “needles”)  targets 
were  used,  Smith  and  Morrison  found  narrow  zones  of  high  uncertainty  for  recovering 
polarizability  for  those  targets  when  the  array  contained  4  single-component  (Z)  receiver  loops. 
These  zones  are  “greatly  ameliorated”  for  the  5-receiver  and  6-receiver  cases  (Fig.  4  in  [3],). 
Although  Smith  et  al  do  study  the  effect  of  receiver  orientation  on  the  recovery  of  target 
polarizability  for  elongated  targets,  they  did  not  study  the  effect  for  vector  (i.e.,  3-C) 
measurements  of  the  secondary  fields. 

4.2. 1.3  Discussion  of  Previous  Receiver  Arrays  Studies 

From  our  review  of  the  two  studies  [2,  4,  5]  we  have  identified  several  important  differences: 

1.  Grimm  [2,  5]  considered  both  single-component  (“1-C”)  and  3-component  (“3-C”) 
receiver  arrays  and  concluded  that  “3-C”  arrays  performed  better.  Smith  et  al  [4] 
considered  only  single-component  receivers. 

2.  Grimm  used  a  set  of  5  different  targets  representing  bodies  with  symmetry  (rod-like  and 
plate-like,  and  sphere-like),  and  two  types  of  asymmetrical  bodies  (clutter-like).  Smith  et 
al  studied  only  two  targets:  a  sphere,  and  a  1000: 1  elongate  target. 

3.  Grimm’s  study  looked  at  the  optimal  number  of  receivers  and  for  this  purpose  he 
evaluated  11  array  geometries  including  the  geometry  of  commercially  available  EMI 
systems  (e.g.  EM-61/EM-63).  Smith  and  Morrison  focused  on  arrays  containing  4,  5,  and 
6  (“1-C”)  receivers  and  studied  performance  of  these  arrays  as  a  function  of  the  receiver 
position.  Although  Grimm’s  study  [2]  did  not  directly  include  the  evaluation  of  an 
antenna  configuration  with  3  lmxlm  transmitter  loops,  he  did  evaluate  a  system 
containing  4  14  x  14  m  transmitter  loops  that  could  be  configured  so  that  the  currents  flow 
with  alternating  polarity  as  well  as  a  common  polarity.  Using  this  configuration  he 
evaluated  performance  based  on  3  different  transmitter  polarities  simulating  in  most 
respects  an  array  with  3  orthogonal  lm  x  lm  transmitters.  The  majority  of  Grimm’s 
testing,  however,  involved  simulating  measurements  over  a  uniform  grid  of  points 
numbering  either  121  (11  x  11  =  Dense  grid)  or  25  (5x5  =  Sparse  Grid).  There  was  no 
effort  in  his  multi-site  modeling  to  introduce  uncertainties  either  in  platform  location  or 
in  platform  attitude.  We  now  know  that  these  uncertainties  dramatically  affect  the  ability 
to  extract  accurate  target  parameters.  Grimm’s  work  concluded  that  the  4-transmitter 
array  performs  significantly  better  than  the  other  10  arrays  he  evaluated  under  all 
scenarios  he  considered. 
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4.2.1.4  About  the  MM  Array 

During  the  AOL,  AOL2,  and  MetalMapper  projects  we  have  made  numerous  free-air  static 
measurements  to  characterize  various  UXO  or  UXO-surrogates  (e.g.,  steel  cylinders).  These 
measurements  were  always  conducted  with  the  targets  centered  beneath  the  antenna  array. 
Furthermore,  the  sites  of  previous  demonstrations  (e.g.,  APG,  and  YPG)  were  flat. 
Consequently,  the  static  data  from  these  sites  was  acquired  with  the  platform  well  centered  over 
the  target.  As  part  of  our  QC  program  at  SLO,  we  acquired  repeat  measurements  over  targets 
anomalies  with  horizontal  offsets  to 
the  source  of  more  than  40  cm. 

Figure  4.9  shows  the  principal 
polarizability  curves  derived  from 
two  MetalMapper  static  data  sets  in 
proximity  to  the  same  target,  a  60 
mm  mortar.  Both  curves  have  good 
Fit  parameters  (99.8  and  99.1).  The 
polarizability  curve  set  from 
SFOStatA01584  (Figure  4.9-  lower 
left)  does  not  exhibit  the  symmetry 
characteristic  we  have  come  to 
recognize  and  even  require  as  an 
essential  attribute  or  indicator  of  a 
UXO-like  shape.  In  this  case,  the 
lateral  offset  between  the  platform 
and  the  target  (~50  cm)  prevented  the 
target  from  being  adequately 
stimulated  along  each  of  its  principal 
axes. 

The  points  we  wish  to  make  here  are: 

1.  Farge  horizontal  offsets  between  a  target  position  and  the  platform  can  occur  in  steep 
terrain.10 

2.  The  symmetry  indicator  derived  from  estimates  of  the  principal  polarizability  can  be 
significantly  degraded  by  large  offsets. 

3.  In  our  opinion,  previous  model  studies  performed  to  select  an  optimum  array  did 
adequately  focus  on  the  effect  of  platform  offset.  Our  intuition  tells  us  that  3 -component 
measurements  of  the  secondary  fields  become  more  important  with  increasing  offset  of 
the  target  from  the  platform. 

4.2. 1.5  Optimal  Arrays  Study  Objectives 

In  the  implementation  of  the  AOF  systems  (AOF  and  AOF2)  as  well  as  the  BUD  system, 
engineering  considerations  led  to  the  construction  of  prototype  antenna  arrays  that  depart  from 
the  optimum  arrays  identified  by  the  cited  studies.  Geometries  is  working  to  commercialize  a 


Yes!  60mm  Mortar 
(Same  target  -  Lateral  Offset) 


Figure  4.9:  An  example  of  the  principal  polarizability  transients 
estimated  from  two  data  sets  over  the  same  UXO  target.  A  large 
lateral  offset  (~0.5m)  between  the  center  of  the  antenna  platform  and 
the  target  in  the  case  of  SLOStatA01584  caused  the  transverse 
polarizability  axes  to  be  inadequately  illuminated.  Note  that  the  shape 
of  the  major  polarizability  transient  (red  curve)  is  nicely  resolved. 


10  Offset  bias  can  be  reduced  if  the  GPS  position  is  corrected  for  platform  attitude  and  heading  in  real-time  before  it 
is  displayed. 
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version  of  this  advanced  EM  technology  based  on  concepts  and  hardware  that  have  evolved 
primarily  out  of  the  AOL  projects  that  were  supported  by  NAVEODTECHDIY  [6].  From  the 
broad  perspective  of  moving  this  technology  from  the  applied  research  and  demonstration  side 
toward  a  point  where  it  is  not  only  accepted  by  the  UXO  and  regulatory  communities  but  also  is 
commercially  available,  it  is  important  that  we  review  the  design  of  the  antenna  arrays  in  order  to 
specify  one  or  more  array  configurations  that  maximize  the  cost-effectiveness  of  deployment  for 
a  particular  purpose. 

Our  objective  here  is  to  expand  on  the  excellent  (and  still  relevant  work)  of  Grimm  and  Sprott 
[2]  and  Smith  and  Morrison  [4]  to  provide  answers  to  the  following  questions: 

1.  Under  what  circumstances  are  full  vector  measurements  of  the  secondary  field  from  a 
target  necessary  or  desirable? 

2.  What  is  the  penalty  that  we  pay  for  reducing  the  number  of  tri-axial  receivers  on  the 
MetalMapper  platform? 

3.  How  does  the  horizontal  offset  between  a  target  and  the  antenna  platform  affect  the 
estimation  of  target  parameters? 

4.  How  sensitive  is  the  performance  of  an  otherwise  optimum  receiver  array  to  small 
changes  in  receiver  positions.  Is  there  a  difference  in  this  sensitivity  for  1-C  and  3-C 
receivers? 

In  this  section,  we  attempt  to  provide  answers  to  these  questions  using  a  new  single-source  target 
parameterization  program,  MMTargets  [1]  and  other  codes. 


4.2.2  Arrays  Studied 

Our  interest  is  in  the  MetalMapper.  Furthermore,  we  are  primarily  interested  in  static 
measurements  at  one  (or  perhaps  more)  sites  that  are  (more-or-less)  over  the  subsurface  target. 
Therefore,  we  use  the  performance  of  the  existing  configuration  of  the  MetalMapper  as  a 
baseline  against  which  we  compare  the 
performance  of  other  systems.  In  Figure  4.10  we 
show  8  arrays  that  we  evaluated.  The  arrays  are 
divided  into  two  groups: 

1.  1-C  -  A  group  of  3  arrays  having  1- 
component  receiver  loops.  In  all  cases, 
these  “1-C”  receiver  loops  are  sensitive 
to  fields  aligned  normal  to  the  plane  of 
the  Z  transmitter. 

2.  3-C  -  A  group  of  5  antenna  arrays 
have  3-component  receiver  cubes. 

Array  names  consist  of  a  4-character  group  plus  a 
2-character  group  separated  by  a  hyphen  (e.g. 

3T5C-30).  The  4  character  group  consists  of  the 
characters  “3T”  that  indicate  the  number  of  transmitter  loops  plus  another  two  characters 
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Figure  4.10:  Antenna  arrays  evaluated  during  this 
study.  All  receiver  coils  are  located  on  the 
diagonals  of  the  Z  (horizontal)  transmitter  coil. 
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indicating  the  number  of  receiver  locations 
and  the  receiver  type  (e.g.,  C  =  Cube  or  tri- 
axial  receiver;  Z  =  single  component 
receiver  sensitive  to  vertical  fields).  The  last 
2-character  group  indicates  the 
MetalMapper  when  it  is  “MM”  and  when 
used  in  connection  with  3T5C  or  3T5Z 
configurations  the  two  figures  indicate  the  X 
and  Y  offset  from  the  center  of  the  loop  to 
the  center  of  the  receiver.  All  receivers  in 
this  study  are  centered  in  the  plane  of  the  Z 
transmitter. 

Note  that,  with  the  exception  of  the 
MetalMapper  configurations  (3T7C- 
MM/3T7Z-MM),  all  the  arrays  have  5 
receivers  (either  C  or  Z)  and  the  only 
difference  between  arrays  is  the  offset  along 
the  diagonal  from  the  array  center.  We  note 
that  the  array  3T5Z-35  approximates  the 
“optimum”  5-receiver  array  determined  by 
Smith  and  Morrison  for  spherical  targets. 


Target  Space 


Primary  ,  Seconary 

Observation  Site  Observation  Site 


Figure  4.11:  Random  target  locations  (250) 
occupying  an  octant  of  a  spheroid  centered  on  the 
horizontal  transmitter  loop  and  having  a  half-length 
equal  to  10X  the  diameter  of  the  target  (100mm  or 
60mm). 


4.2.3  Target  Space 

Given  the  constraints  on  the  array  geometries  (i.e.,  receivers  lie  along  the  diagonals  of  the 
transmitter  loops),  we  choose  to  constrain  all  targets  to  an  octant  of  a  spheroid  that  is  centered  at 
the  center  of  the  Z  transmitter  loop  (also  the  center  receiver)  having  a  half-length  equal  to  the 
10X  the  diameter  of  the  test  target  (e.g.,  lm  for  a  100mm  sphere)  and  a  radius  equal  to  half  the 
side  dimension  of  the  Z  transmitter  loop  (0.5m).  Plan  and  cross-sectionals  of  random  target 
locations  appropriate  for  each  target  diameter  are  shown  on  the  left  (100mm  Sphere  target)  and 
on  the  right  (60mm  Prolate  Spheroid)  of  Figure  4.1 1.  The  plotted  points  represent  the  locations 
of  250  random  targets.  Our  experience  with  the  MetalMapper  based  on  demonstrations  and  a 
limited  amount  of  modeling  is  that  parameter  extraction  is  much  more  reliable  when  the  targets 
have  a  horizontal  offset  that  is  less  than  0.4m  ( r<0.4m ).  In  Figure  4.1 1,  we  have  drawn  a  second 
(inner)  boundary  that  corresponds  to  this  smaller  (0.4m)  offset  distance.  We  call  the  targets  that 
are  confined  within  this  smaller  spheroid  octant  the  Reduced  Target  Set  (RS). 

Our  single-source  dipole-based  modeling  program  (MMTargets)  is  able  to  extract  target 
parameters  from  data  measured  at  one  or  more  measurement  sites.  In  this  report,  we  include  a 
few  results  based  on  the  analysis  of  data  from  2  sites.  The  primary  site  (shown  with  a  magenta 
disk)  is  at  the  center  of  the  plan  view  of  target  space  (Figure  4.1 1).  The  secondary  site  is  offset 
by  0.5m  east  and  therefore  is  centered  and  tangent  to  the  east  edge  of  the  full  target  set  (FS).  We 
have  chosen  to  restrict  our  target  placement  to  a  single  octant  of  a  spheroidal  volume.  Using 
solutions  based  on  a  measurement  at  the  primary  measurement  site,  symmetry  arguments  allow 
us  to  infer  that  the  results  would  be  the  same  if  we  distributed  the  targets  over  the  other  3  octants 
and  processed.  However,  to  make  the  same  inferences  using  a  multi-site  solver  we  would  need 
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to  occupy  5  sites  (a  center  site  plus  4  sites  offset  in  the  4  cardinal  directions).  Here,  we  are 
simulating  a  case  where  the  location  of  the  target  based  on  a  measurement  at  the  primary  site 
vectors  us  to  the  secondary  site. 


4.2.4  Target  Features 

In  this  study  we  follow  the  lead  of  Smith  and  Morrison  [4]  and  use  a  100mm  sphere  as  the  target 
for  the  initial  evaluation  of  array  performance.  Then  we  check  our  results  with  a  target  that 
better  simulates  the  response  of  a  UXO.  For  this  purpose,  we  use  a  60mm  diameter  prolate 
spheroid  shape  with  a  3:1  aspect  ratio  (length  =  180mm).  The  principal  polarizability  transients 
of  the  two  targets  are  plotted  in  Figure  4.12.  For  the  purpose  of  this  study,  the  polarizability  is 
sampled  at  42  time  logarithmically  spaced  time  gates  ranging  from  106<t<7916  ps.  Each  time 
data 


Principal  Polarizability  Curves 
Array  Study  —  Target  Objects 

Principal  Polarizability  Curves 
lOOirnn  Sphere 


100  200  500  1000  2000  5000  lxlO4 

Time  Ips) 

Principal  Polarizability  Curves 
60inin  Prolate  Spheroid  (3:1  Aspect  Ratio) 


Table  4.1:  Tabulation  of  various 
target  features  extracted  from 
principal  polarizabilities  shown  in 
Figure  6  and  the  target  dimensions. 


we  invert  a  data  set, 

MMTargets  provides  an 
estimate  of  the  principal 
polarizability  transients  of 
the  target  together  with  its 
extrinsic  properties  (i.e., 
the  position  and  attitude 
angles).  In  principle,  one 
could  match  the  resulting 
estimated  target 

polarizabities  with  the 
appropriate  curves  in 
Figure  4.12  and  form  an 
error  estimate  (e.g.,  MSE, 

RMSE,  etc).  In  our  case,  it 
was  easier  to  extract  a  few 
scalar  features  from  the 
curves  and  compare  those 
with  known  values 

determined  either  theoretically  or  calculated  directly  from 
the  curves  in  Figure  4.12.  We  chose  the  latter.  In  Figure 
4.13  we  define  a  set  of  5  meta  features  that  we  extract  from  the  parameters  that  are  provided  by 

the  modeling  program  (MMTargets).  One  feature  (V)  is  simply  proportional  to  the  spheroid 
volume  (estimated  from  the  spheroid  length  and  diameter).  The  other  4  features  are  calculated 
from  the  integration  of  the  0th  and  1st  order  time  moments  of  the  3  principal  polarizability  curves. 
The  true  values  of  these  parameters  for  the  100mm  sphere  and  the  60mm  prolate  spheroid  target 
have  been  tabulated  in  Table  4. 1 . 


100  200  500  1000  2000  5000  lxlO4 

Time  (fa) 

Figure  4.12:  Principal  polarizability 
transients  for  the  100mm  sphere 
(top)  and  the  60mm  diameter  prolate 
spheroid  used  as  targets  in  this 
study. 
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Figure  4.13:Definition  of  5  target  features  used  in  this  report  as  the  basis  for  an  index  of 
array  performance.  The  values  of  these  5  parameters  are  known  (Table  2.1).  The 
distribution  of  the  parameter  estimates  in  linear  or  log  space  suggests  that  they  form  a 
Gaussian  normal  distribution. 


4.2.5  Electromagnetic  Interference 

The  array  studies  referenced  earlier  [2, 

4,  5]  corrupted  otherwise  perfect  model 
data  with  random  noise  to  make  the 
array  study  more  realistic.  Grimm  [2] 
based  the  amplitude  of  his  noise  levels 
on  a  fraction  of  the  signal  level  of  his 
cylinder  target  oriented  vertically  and 
placed  at  the  shallowest  depth 
(0.55m11).  Smith  [4],  on  the  other 
hand,  assumed  a  constant  noise  field  of 
2  nT/s  for  receivers  responding  to  a 
vertical  field,  varying  to  6  nT/s  for 
receivers  responding  to  horizontal 
fields.  Our  noise  data  suggest  that 
their  estimate  of  noise  is  low  by  as  much  as  an  order  of  magnitude  depending  on  stacking.  In 
this  study,  we  base  our  noise  values  on  noise  measurements  that  we  have  derived  from 
experimental  measurements  with  the  MetalMapper  over  background  points.  We  show  in  Figure 


I  Time  (ms)  Time  (ms)  I 

Figure  4.14:  :  RMS  noise  levels  measured  with  the  MM  at  Camp 
Butner,  NC  in  2010.  The  noise  levels  are  shown  for  the  case 
when  the  Z  transmitter  is  energized.  The  noise  when  the  X  and 
Y  transmitters  are  energized  is  similar  in  magnitude  and 
behavior. 


11  The  depth  quoted  here  includes  0.4m  which  is  the  assumed  height  of  the  horizontal  transmitter  loop  above  ground- 
level. 
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4.14  transients  representing  the  RMS  noise  values  for  each  time  channel  and  each  receiver  the 
case  of  no  stacking.  The  results  are  plotted  only  for  the  case  where  the  Z  transmitter  has  been 
energized.  Plots  made  with  the  X  and  Y  transmitter  energized,  respectively,  produce  plots  that 
have  RMS  noise  values  similar  both  in  magnitude  and  decay  rate  with  those  shown  for  the  Z 
transmitter  (Figure  4.14).  The  decay  in  noise  values  in  time  is  characteristic  of  the  behavior  of 
RMS  noise  levels  when  presented  this  way  since  it  reflects  the  greater  averaging  that  arises 
because  of  the  logarithmically  increasing  gate  widths  as  the  gate  centers  are  moved  out  in  time. 
Figure  4.14  suggests  that  the  noise  levels  at  delay  times  of  around  0.1  ms  (lOOus)  are  on  the 
order  of  100  nT/s.  The  noise  level  is  larger  by  about  a  factor  of  2  for  the  horizontal  components. 
Therefore,  we  have  added  random  noise  to  our  synthetic  data  in  the  amount  of  300  nT/s  for  both 
horizontal  components  and  150  nT/s  on  the  vertical  components. 12  Consistent  with  the 
logarithmically  expanding  gate  widths  at  later  times,  we  reduce  the  level  of  noise  at  each 
subsequent  gate  by  a  factor  equal  to  the  square  root  of  the  ratio  of  the  gate  width  at  106us  to  the 
gate  width  of  the  time  gate  under  consideration.  In  practice,  this  causes  the  noise  level  to 
decrease  with  a  slope  of  -  1/2  when  plotted  in  Log-Log  space  as  is  done  in  Figure  4.14. 

4.2.6  Position  and  Attitude  Uncertainty 

When  processing  static  data  from  a  single  measurement  site,  small  errors  in  platform  position 
and  attitude  angles  simply  serve  to  bias  the  estimated  position  and  attitude  of  the  target  when 
they  are  transformed  into  the  geographic  coordinate  system.  However,  we  computed  some 
results  using  two  measurement  sites.  In  order  to  show  how  small  errors  in  position  and  platform 
attitude  affect  multi-site  performance,  we  introduced  random  variations  in  each  of  the  three 
coordinates  with  standard  deviation  of  0.02m  and  in  the  3  attitude  angles  with  standard  deviation 
1°  when  calculating  the  model  data. 

4.2.7  Parameter  Extraction 

We  inverted  each  set  of  250  data  files  using  MMTargets,  the  single-source  modeling  program 
that  has  been  described  in  some  detail  in  Section  2  of  this  report.  As  related  in  that  section, 
MMTargets  generates  constrained  solutions  for  the  best  fitting  sphere,  prolate  spheroid,  and 
oblate  spheroid.  Lastly,  it  solves  the  fully  unconstrained  problem  for  the  best  fitting  ellipsoid. 
From  those  solutions,  we  extracted  the  position  and  attitude  estimates  together  with  the  estimates 
for  the  target’s  radius  and  half-length.  In  a  separate  operation,  MMTargets  computes  the 
principal  polarizability  transients  using  the  target  position  corresponding  to  the  constrained 
solution  that  has  the  least  mean  square  error  (MSE).  Predictably,  that  solution  invariably 
corresponds  to  the  ellipsoidal  model.  From  these  target  parameters  (primary  features)  we  are 
able  to  compute  the  5  meta-parameters  (and  more  that  are  not  listed)  highlighted  in  Figure  4.13. 
In  subsequent  sections  we  will  use  scatter  plots  of  scalar  parameters  that  we  have  compiled  from 
each  the  data  sets  corresponding  to  each  array  and/or  target  situation.  And  finally,  we  will  show 
a  metric  that  quantifies  in  some  sense  the  relative  performance  of  various  arrays. 


12  In  retrospect,  we  realized  that  the  noise  levels  shown  in  Figure  8  do  not  account  for  stacking.  Therefore,  we  could 
justify  reducing  the  noise  levels  by  a  factor  of  10  to  account  for  the  effect  of  stacking.  This  would  no  doubt  greatly 
improve  our  performance  estimates.  But  the  noise  levels  would  remain  at  10/20  nT/s,  which  is  several  times  larger 
than  that  used  by  Smith  [3], 
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4.2.8  SNR  and  MSE 

In  general,  the  quality  of  the  target 
parameters  derived  from  a  data  set 
corresponding  to  a  static  MetalMapper 
measurement  of  a  target  located  at  a 
specified  position  with  a  specified  attitude 
can  be  determined  by  the  quality  of  the  fit 
to  the  observed  data.  In  MMTargets,  the 
fit  parameter  is  the  mean  square  error 
(MSE)  between  the  input  data  and  the  data 
calculated  based  on  estimates  of  the 
various  target  parameters.  In  Figure  4.15, 
we  summarize  data  quality  with  two 
histograms  and  two  scatter  plots  using  the 
SNR  and  the  MSE  fit  statistics  for  the  set 
of  250  sphere  data  sets  corresponding  to 
the  3T7C-MM  (the  MetalMapper  array) 
and  the  relationship  between  the  two  quantities.  A  few  remarks  are  in  order.  Firstly,  the 
histogram  of  SNR  shows  that  about  80  percent  of  the  data  sets  processed  have  an  SNR  >  20dB. 
Secondly,  these  uniformly  high  SNR’s  result  in  a  similar  percentage  of  the  corresponding 
MMTargets  solutions  exhibiting  rather  good  MSE  (i.e.,  MSE<0.2)  as  shown  by  the  histogram  of 
the  MSE’s  for  this  data  set  (Figure  4.15B).  In  the  right  hand  side  of  Figure  4.15(panels  C  &  D) 
we  show  scatter  plots  for  the  full  (250  targets)  target  data  set  (Figure  4.15C)  and  the  reduced 
target  set  (164  targets).  Recall  that  in  the  reduced  set  (RS)  we  have  excluded  targets  positioned 
outside  of  the  inner  spheroidal  boundary  (see  Figure  4.1 1).  What  is  interesting  about  these  plots 
is  that  the  outlier  points  shown  in  Figure  4.15C  (FS)  disappear  when  we  reduce  the  target  set  by 
removing  targets  outside  of  the  spheroid  of  radius  0.4m  (see  Figure  4.11).  This  implies  that 
MMTarget  has  trouble  finding  good  solutions  for  targets  with  horizontal  offsets  that  place  them 
beneath  the  sides  of  the  transmitter.  Detailed  examination  of  shallow  targets  with  horizontal 
offsets  in  the  range  0.4<r<0.5m  show  many  if  not  all  the  targets  forming  the  outlier  points  in 
Figure  4.15C  fall  within  this  zone. 

4.2.9  Target  Size  &  Time  Persistence 

Target  parameters  related  to  size  and  to  time  persistence  are  frequently  very  useful  target 
parameters  for  classification.  In  Figure  4.13,  we  have  identified  2  size-related  target  features 

( V ,  ho)  and  a  time  persistence  feature  {zn)  that  we  will  use  to  evaluate  array  performance.  In 
Figure  4.16,  we  show  Log-Log  scatter  plots  of  these  parameter  pairs  for  the  full  set  (FS)  of 
targets  and  for  the  reduced  set  (RS).  These  scatter  plots  show  that  both  size  estimates  and  the 
time  estimates  cluster  well  around  their  expected  values  in  Log  space.  We  have  also  included 
histograms  showing  the  distribution  of  the  zn  (in  linear  space  -  bottom  left  of  Figure  4.16)  and 
the  two  size  features  plotted  in  logarithmic  space.  To  emphasize  that  the  mean  of  these 
parameters  effectively  equal  to  their  theoretical  values,  the  appropriate  value  has  either  been 
subtracted  (as  in  the  case  for  zn  -  Figure  4.16  lower  left)  or  used  as  a  normalization  factor  (the 
two  size  parameters  -  bottom  center  and  right  panels  in  Figure  4.16).  These  histograms  show  that 
time  errors  approximate  a  distribution  with  (approximate)  bilateral  symmetry  in  linear  space 
while  the  two  size  parameters  approximate  distributions  with  bilateral  symmetry  in  log  space 


B)  Distribution  of  MSE  D)  Reduced  Target  Set  (RS) 

Figure  4.15:Summary  of  the  distribution  of  SNR  (A),  MSE 
(B)  and  the  relationship  between  the  two  (C  &  D).  The 
figures  are  based  on  inversions  of  250  target  data  sets  for  the 
3T7C-MM  array. 
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Figure  4.16:  Scatter  plots  showing  target  clustering  based  on  2  different  size  parameters  (Size  A/SizeB).  For  each  of 
these  parameters  we  have  used  the  time  persistence  (tI2)  as  the  plot  ordinate.  The  histograms  shown  below  the  scatter 
plots  exhibit  approximate  bilateral  symmetry  about  a  mean  value  of  0.  Thus  a  good  measure  of  array  performance 
based  on  these  3  parameters  would  be  the  standard  deviation  (o)  of  the  plotted  parameter. 

(e.g.,  log-normal  distributions).  We  emphasize  this  behavior  since  we  will  use  the  standarc 
deviations  of  these  parameters  as  performance  metrics  for  the  arrays  we  evaluate. 


4.2.10  Target  Shape 

After  size  and  time  persistence,  target  attributes  that  indicate  shape  are  particularly  important.  In 
fact,  many  times  interpreters  (this  author  included)  focus  on  target  symmetry  as  indicated  in  a  set 
of  principal  polarizability  curves  having  a  single  major  polarizability  transient  and  a  pair  of 
smaller  minor  transients  that  are  identical  or  nearly  identical.  This  is  the  ideal  case  and  the  one 
that  we  always  see  in  free-air  or  pit  measurements  when  the  target  is  located  beneath  the  center 
of  the  MetalMapper  antenna  array.  But  as  Figure  4.9  illustrates,  one  should  be  careful  to  take 
into  account  the  offset  of  the  target  from  the  center  of  the  array.  We  now  know  that  large 
horizontal  offsets  can  degrade  target  features  affected  by  the  targets  symmetry.  In  Figure  4.9,  the 
lack  of  symmetry  suggested  by  the  polarizability  curve  set  shown  in  the  panel  on  the  lower  left  is 
caused  by  the  fact  that  the  target  is  offset  from  the  center  of  the  MetalMapper  array  by  a  distance 
of  slightly  more  than  0.5m. 
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In  Figure  4.13  we  have  identified  2  target  attributes  that  provide  important  information  about 
target  shape: 

1.  Por:  This  is  a  measure  of  the  aspect  ratio.  We  have  defined  it  as  the  ratio  of  a  scalar 
measure  of  the  major  polarizability  curve  to  the  geometric  mean  of  the  corresponding 
scalar  measures  of  the  two  minor  polarizability  curves.  Note  that  in  our  definition, 
aspect  ratio  is  independent  of  symmetry.  We  expect  a  sphere  to  have  an  aspect  ratio  of 
1.  A  UXO  is  generally  expected  to  be  elongated  and  therefore  will  have  an  aspect 
ratio  P or>  1 . 

2.  RLo:  We  use  this  parameter  as  an  indicator  of  target  symmetry.  The  parameter  is 
formed  as  the  ratio  of  two  size-related  scalar  metaparameters  derived,  respectively, 
from  the  two  minor  principal  polarizability  transients.  When  the  minor  curves  are 
identical,  RLo=l.  Because  we  order  the  principal  polarizability  according  to  their  size 
on  some  basis,  good  data  will  usually  produce  values  that  are  greater  than  or  equal  to  1 
(RL0>1).13 


Using  the  parameters  extracted  from  our  sphere  data  (3T7C-MM  array),  we  have  studied  these 
two  shape  parameters  and  their  dependence  on  distance  (R)  from  the  center  of  the  MetalMapper 
array,  the  horizontal  offset  distance  (r)  from  the  array  center,  and  the  target  depth.  Scatter  plots 
for  those  parameter  pairs  are  in  Figure  4.17A  and  Figure  4.17B. 

In  the  top  row  of  Figure  4.17A,  the  aspect  ratio  is  plotted  against  target  distance  or  range  (R). 
When  looking  at  the  full  set,  we  see  that  estimated  aspect  ratios  show  significant  departure  from 
their  true  value  (1)  at  a  range  of  0.5  and  then  again  at  ranges  near  the  maximum  for  those  points 
plotted.14  Note  that  the  outlier  points  at  R=0.5m  are  gone  when  we  plot  the  points  for  the 
reduced  set.  One  can  infer  from  these  2  plots  alone  that  the  outliers  come  from  shallow  targets 
with  horizontal  offsets  near  0.5m.  The  second  and  3rd  pair  of  plots  confirm  that  the  problem 
targets  are  shallow  and  have  offsets  near  0.5m. 

In  Figure  4.17B,  we  use  RLo  as  the  shape 
parameter.  The  plot  pairs  once  again  show 
that  our  ability  to  determine  symmetry 
properties  of  targets  degrades  both  with 
horizontal  offset  (middle  pair)  and  with 
depth  (lower  pair).  We  are  clearly  attuned 
to  the  fact  that  the  quality  of  the  target 
parameters  degrades  with  low  SNR  and  so 
the  scatter  in  the  symmetry  plots  at  greater 
target  depths  comes  as  no  surprise.  For 


13  Our  definition  of  RL0  is  based  on  numerical  integrations  of  the  principal  polarizability  transients.  There  are  cases 
usually  involving  noisy  data  or  polarization  curves  exhibiting  cross-over.  In  these  cases,  we  sometimes  see  RLo<l. 

14  The  astute  reader  will  notice  that  the  range  (R)  of  the  target  sets  halts  at  R=0.8m  and  not  1.0m.  That  is  because 
we  used  a  noise  level  that  was  too  large  to  get  good  solutions  for  the  deeper  targets.  We  excluded  those  from  these 
plots. 
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completeness  and  in  order  to  justify  the  use  of  these  two  shape  parameters  in  defining  a 
performance  metric,  we  include  histograms  showing  the  distribution  of  the  both  Por  and  RLo  as 
Figure  4.18. 


4.2.11  An  Array  Performance  Metric 

The  preceding  discussions  together  with  the  various  graphics  have  been  leading  up  to  a  definition 
of  a  metric  function  that  we  can  use  to  summarize  the  performance  of  the  arrays  that  we  evaluate. 
That  end,  the  five  parameter  histograms  in  Figure  4.16  (bottom  panel)  and  Figure  4.18  suggest 
that  the  standard  deviations  of  the  parameters  plotted  can  be  used  as  a  metric.  In  equation  4.9 
below,  we  define  that  metric. 


(  )  =  Standard  Deviation 


(4.9) 


In  equation  4.9,  the  5  quantities  ( V,I20,POE,RLo,r1 )  are  the  estimates  of  the  5  features  for  the 
target  population  for  a  given  array  situation.  The  two 
scalar  parameters  (V,I2)  are  the  known  values  of  the 

spheroid  volume  and  the  scalar  polarizability  tensor 
invariant  for  the  target  being  used  (see  Table  4.1).  We 
have  chosen  to  compute  standard  deviations  of  3  target 
features  after  taking  their  logarithms  because  these 
parameters  appear  bilaterally  symmetric  and  have 
mean  values  at  or  near  their  correct  values. 

4.2.12  Array  Performance  Results 

As  we  indicated  earlier  in  this  report,  we  first  evaluated 
the  performance  of  6  of  the  arrays  (the  5  3-C  arrays 
plus  3T7Z-MM)  shown  in  Figure  4.10.  The  results 
reflect  inversions  from  a  single  receiver  site.  Then, 
using  the  multiple-site  inversion  capability  of 
MMTargets,  we  studied  the  performance  of  the  present 
MetalMapper  array  when  two  sites  over  the  target  have 
been  measured.  In  this  case,  we  simulate  the  scenario 
of  a  repeated  measurement  based  on  large  target  offset. 

In  such  a  case,  the  repeat  measurement  site  has  been 
offset  to  a  site  based  on  the  estimated  target  position 
using  data  from  an  earlier  measurement  site. 


Performance  —  6  Array  Candidates 
Sphere  Target  (FS) 


Performance  —  6  Array  Candidates 
Sphere  Target  (RS) 


Figure  4.19:  Performance  summary  for  6 
arrays  evaluated  using  the  sphere  target  set. 
The  upper  panel  represents  results  for  the  full 
set  (FS),  the  lower  panel  shows  results  for  the 
reduced  (RS)  target  set. 


4.2.13  Single-Site  Data  -  Sphere  Target 

The  performance  results  for  all  5  arrays  containing  tri-axial  receivers  (3-C)  plus  one  array 
containing  only  single-component  (1-C)  receivers  are  summarized  by  means  of  stacked  bar 
graphs  in  Figure  4.19.  Perfect  performance  would  result  if  the  standard  deviations  computed  in 
equation  4.9  were  zero.  Therefore  performance  increases  in  the  reverse  direction  of  the  P  value. 
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For  both  the  full  and  reduced  data  sets,  the  current  MetalMapper  configuration  (3T7C-MM) 
performs  the  best.  Next,  notice  that  the  1-C  version  of 
the  MetalMapper  configuration  (3T7Z-MM)  is  only 
marginally  better  than  the  array  with  the  worst 
performance  score  (3T5C-20).  Clearly,  3-component 
receivers  significantly  improve  performance,  even  when 
the  targets  are  kept  well  within  the  lm  x  lm  area 
defined  by  of  the  horizontal  transmitter.  Finally,  there 
is  a  clear  trend  that  suggests  that  in  a  receiver  array 
consisting  of  5  Cubes,  the  receivers  should  be  placed 
outside  the  geometric  center  of  each  quadrant  of  the  lm 
x  lm  transmitter  loop.  The  results  are  slightly 
ambiguous  on  that  score  since  the  3T5C-35  array  is  the 
best  5 -Cube  array  when  evaluated  with  the  full  set  (FS) 
of  sphere  targets.  But  when  evaluated  against  the 
reduced  set  (RS),  the  3T5C-30  performs  better.  We 
note  here  that  in  their  study  Smith  et  al  [4]  determined 
that  a  system  with  5  1-C  receivers  constrained  to  be 
placed  within  the  transmitter  loop  is  optimum  for  sphere 
targets  when  the  receivers  are  near  the  center  and  near 
the  4  comers  of  the  transmitter  loop  at  an  offset  distance 
of  between  50  and  60  cm.  The  3T5C-35  has  4  of  its  3-C 
receivers  at  a  distance  of  -50cm  from  the  center.  The 
3T7C-MM  array  has  2  receivers  at  a  distance  of  55  cm. 

Our  results  appear  to  be  consistent  therefore  with  the  findings  of  Smith  and  his  associates. 

4.2.14  Single-Site  Data  -  Prolate  Spheroid  Target 

Using  the  60mm  prolate  spheroid  (PS)  as  a  target  model,  we  evaluated  the  performance  of  the  3 
arrays  showing  the  best  performance  based  on  the  sphere  target  (3T5C-30,  3T5C-35,  and  3T7C- 
MM).  We  present  those  results  in  Figure  4.20.  As  with  the  sphere  target,  we  have  evaluated  the 
arrays  for  both  the  full  target  set  (Figure  4.20  top  panel)  and  for  the  reduced  target  set  (Figure 
4.20  bottom  panel). 

The  results  in  Figure  4.20  confirm  with  a  second  target  geometry  that  the  current  MetalMapper 
array  geometry  is  better  than  the  best  of  the  arrays  having  5  3-C  receivers.  Note,  however,  that 
on  the  reduced  target  set  (Figure  4.20  bottom  panel)  the  overall  performance  of  the  3  arrays  is 
almost  identical  with  the  MetalMapper  maintaining  a  slight  edge.  We  think  this  is  due  to 
improved  performance  of  the  MetalMapper  on  those  targets  having  larger  horizontal  offsets. 
And  notice  too  that  the  difference  in  performance  comes  largely  from  the  ability  to  extract  better 
estimates  of  the  “Shape  A”  parameter  (target  aspect  ratio  Por)  with  data  from  the  MetalMapper 
array. 

We  have  combined  the  performance  of  the  3  arrays  for  the  sphere  target  set  (Figure  4.19  and  the 
prolate  spheroid  target  set  (Figure  4.20)  into  a  set  of  bar  graphs  summarizing  the  performance 
when  considering  both  targets.  Those  results  are  shown  in  Figure  4.21.  Here  again,  the  3T7C- 
MM  clearly  retains  its  performance  advantage  against  the  2  competing  5  cube  systems, 
particularly  when  considering  the  full  target  sets.  But  when  using  the  reduced  target  set  and  thus 
eliminating  all  targets  with  offsets  outside  the  0.4m  limit,  the  5-Cube  arrays  do  better  with  the 


Performance  —  3  Array  Candidates 
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Performance  —  3  Array  Candidates 
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Figure  4.20:  Performance  summary  for  the  3 
best  3-C  arrays  from  evaluation  with  the 
sphere  target.  The  performance  is  based  on 
parameters  extracted  from  the  60mm  prolate 
spheroid  target  set  with  single-site  parameter 
extraction.  The  top  panel  represents 
performance  with  the  full  target  set.  The 
bottom  panel  is  for  the  reduced  target  set. 
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3T5C-30  array  performing  marginally  better  than  the  3T5C-35  system  that  has  its  receivers  7  cm 
farther  away  from  the  center. 


4.2.15  Multi-Site  Performance 

In  principle,  measurements  taken  at  several  positions  (“sites”)  over  the  targets  should  vastly 
improve  the  fidelity  (i.e.,  performance)  of  the  modeling  program  used  target  parameter  extraction 
when  the  program  is  capable  of  processing  multi-site  data  sets.  MMTargets  has  multi-site 
inverse  modeling  capability  for  a  single  (dipole)  target  source.  We  used  the  multi-site  capability 
to  study  performance  under  the  scenario  that  we  have  experienced  repeatedly  in  our 
demonstrations  at  San  Luis  Obispo,  Butner,  and  most  recently  at  Mare  Island.  The  scenario  is 
that  our  first  measurement  indicated  that  the  target  is  offset  by  more  than  0.4m  from  the  center  of 
the  MetalMapper  array.  A  second  data  site  was  placed  at  the  target  position  estimated  from  the 
original  measurement.  Thus,  we  have  data  available  over  (hopefully)  a  single  target  at  2 
different  surface  sites. 


We  evaluated  the  performance  of  2-Site  data  sets  using  the  sphere  target  and  the  MetalMapper 
array.  The  two  target  sites  are  noted  in  Figure  4. 1 1  (“Primary”  and  “Secondary”  observation 
sites).  The  primary  site  is  located  at  the  center  of  the  spheroid  whose  octant  volume  contains  the 
random  target  positions.  The  secondary  site  is  offset  laterally  along  the  X  axis  by  a  distance  of 
0.5m.  Together  the  measurements  at  the  two  site  assure  that  one  of  the  sites  will  be  reasonably 
centered  over  the  target.  We  evaluated  performance  for  2  cases: 

1 .  Perfect  Position  and  Attitude  -  In  this  case,  the  position  and  antenna  attitude  angles 
are  assumed  to  be  known  perfectly. 

2.  Reasonable  Position  and  Noise  Uncertainty  -  In  this  case,  we  allowed  for  random 
errors  in  both  the  platform  position  (a=2cm  for  all  3  coordinates)  and  random  errors  in 
the  3  attitude  angles  (o=l°).  The  standard  deviations  used  for  these  uncertainties  are 
consistent  with  those  reported  for  RTK  GPS  positioning  and  for  the  module  that 
transduces  attitude  angles  for  the  sensor  platform.15 
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Figure  4.21:  Performance  summary  for  the  3  best  3-C  arrays  from  evaluation  with  the  sphere  target.  The 
performance  is  based  on  parameters  extracted  from  both  the  100mm  sphere  and  the  60mm  prolate  spheroid  target 
sets  with  single-site  parameter  extraction.  The  left  panel  represents  performance  with  the  full  target  set  (FS).  The 
right  panel  is  for  the  reduced  target  set  (RS). 


15  Platform  pitch  and  roll  angles  have  specified  accuracies  on  the  order  of  ±0.1°.  Unfortunately,  the  measurement 
error  for  the  heading,  which  are  based  on  vector  measurements  of  the  magnetic  field,  are  an  order  of  magnitude 
larger.  We  think  that  random  angle  errors  having  a  standard  deviation  o=l°  may  in  fact  be  optimistic. 
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The  performance  results  for  2-site  parameter  extractions  are  shown  in  Figure  4.22  as  stacked  bar 
graphs.  These  results  are  shown  for  the  3T7C-MM  and  the  3T5C-30  arrays.  For  comparison,  we 
have  also  included  the  performance  (middle  bar)  for  single-site  parameter  extraction  for  those 
arrays.  These  results  show  that  2-Pt  and,  more  generally,  multi-point  parameter  extraction  can 
significantly  improve  performance.  But  this  improvement  can  only  be  achieved  when  the 
precision  of  both  platform  attitude  angles  and  positioning  are  perhaps  an  order  of  magnitude 
smaller  than  the  uncertainties  that  we  used  for  these  evaluations. 
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Figure  4.22:  Performance  results  for  2-point  parameter  extraction  with  the  3T7C-MM  (Left  panel)  and  3T5C-30 
(Right  panel)  arrays.  The  performance  estimates  uses  parameter  estimates  from  the  sphere  target  group.  For 
comparison,  we  show  the  performance  of  the  single-point  extractions  (center  bar). 


4.2.16  Conclusions-  Optimum  Array  Study 

In  this  study  we  have  evaluated  the  performance  of  the  existing  MetalMapper  array  (357C-MM) 
against  a  set  of  candidate  arrays  having  fewer  tri-axial  receivers  (5).  We  have  also  investigated 
whether  there  is  a  significant  benefit  in  deploying  tri-axial  receivers  rather  than  a  set  of  vertical 
receiver  loops  such  as  used  in  the  BUD  array.  We  emphasize  that  our  conclusions  relate  to  the 
efficacy  of  the  arrays  for  static  measurements.  However,  earlier  studies  that  we  have  cited 
strongly  suggest  our  conclusions  extrapolate  to  the  case  of  dynamic  measurements  with  an  array 
containing  a  single  Z  transmitter  and  an  array  of  tri-axial  receivers.  Our  results  lead  to  the 
following  conclusions: 

1.  The  MetalMapper  with  its  present  7-Cube  receiver  array  performs  significantly  better 
than  the  5 -Cube  symmetrical  receiver  arrays  regardless  of  receiver  spacing.  However, 
if  we  restrict  our  target  space  to  a  volume  such  that  shallow  targets  are  never  offset 
horizontally  farther  than  0.4m  from  the  array  center,  the  array  3T5C-30  has 
performance  that  is  not  significantly  degraded  from  that  of  3T7C-MM.  We  believe 
that  the  two  tri-axial  cube  receivers  located  at  coordinates  (0.39,0.39)  and  (-0.39,-0.39) 
or  about  55cm  from  the  center  of  the  array  account  for  the  improved  performance  of 
the  MetalMapper  array  over  its  5-Cube  rivals. 

2.  The  results  show  clearly  that  if  we  restrict  ourselves  to  5  tri-axial  receivers,  we  should 
place  those  receivers  along  the  diagonals  at  a  distance  of  between  42cm-55cm  from 
the  array  center.  Our  results  show  that  3T5C-30  (42  cm  diagonal  distance  from  center) 
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performs  marginally  better  3T5C-35  (49.5cm  diagonal  distance  from  center).  Such  a 
wide  spacing  is  consistent  with  the  optimum  1-C  receiver  placement  as  determined  by 
Smith  et  al  [4], 

3.  We  conclude  that  there  is  a  substantial  benefit  to  making  tri-axial  measurements.  This 
conclusion  is  made  on  the  basis  of  the  poor  performance  of  the  3T7Z-MM  array 
relative  to  its  3-C  counterpart  (3T7C-MM).  Although  we  were  not  able  to  study  the 
performance  of  the  2  other  candidate  1-C  arrays,  we  have  no  reason  to  believe  that 
fewer  1-C  receivers  would  surprise  us  by  performing  better  than  their  3-C  counterparts. 
We  note  that  array  3T7Z-MM  performs  better  on  the  reduced  target  set.  So  in  all 
probability,  if  we  wish  to  use  1-C  receivers,  the  penalty  we  will  pay  is  that  we  will 
need  to  reduce  the  maximum  permissible  horizontal  offset  of  the  target  from  the 
antenna  platform. 

4.  Our  results  show  that  the  performance  of  these  arrays  is  generally  well  characterized 
by  using  spherical  targets.  But  there  is  some  ambiguity  in  those  results  because  the 
ranking  of  3T5C-30  and  3T5C-35  depend  on  whether  we  are  using  the  full  set  or  the 
reduced  set.  However,  when  we  composite  the  results  from  both  the  sphere  target  and 
the  prolate  spheroid  the  ranking  of  the  3  arrays  is  consistent  between  the  full  set  and 
the  smaller  reduced  set. 

5.  Admittedly,  our  efforts  to  evaluate  the  efficacy  of  multi-site  parameter  extraction  were 
limited.  Nevertheless,  our  experiments  suggest  that  the  performance  of  multi-site 
parameter  extraction  will  be  inferior  to  that  of  a  single  well  sited  measurement  over 
the  target  unless  we  develop  an  efficient  means  for  locating  the  platform  with  sub¬ 
centimeter  accuracy  and  transducing  its  3  attitude  angles  with  0.1°  accuracy.  This 
level  of  (internal)  accuracy  might  be  achieved  by  either  moving  the  antenna  array  to 
precisely  located  positions  on  a  rigid  planar  template,  or,  alternatively,  employing  an 
array  (e.g.  2  x  2)  of  identical  antenna  arrays.  Such  an  array  would  be  similar  in 
concept  to  the  TEMTADs  array. 


4.2.17  Recommendations  for  Further  Study 

There  are  a  number  of  areas  where  further  study  would  be  useful. 

1.  Our  target  sets  were  generated  with  a  single  noise  level  consistent  with  noise 
measurements  at  several  demonstration  sites.  During  analysis  we  discovered  that  we 
did  not  account  properly  for  the  noise  reduction  due  to  stack  measurements  (VN).  As  a 
result,  the  noise  levels  were  perhaps  overly  high.  The  high  noise  does  not  affect  the 
relative  rankings  of  the  arrays  since  data  sets  for  each  array  were  generated  with 
identical  noise  statistics.  However,  we  would  like  to  investigate  different  noise  levels 
and  how  they  affect  the  overall  results. 

2.  We  would  like  to  extend  the  multiple  site  study  in  2  ways.  We  would  like  to 
experimentally  determine  the  required  specification  for  both  position  and  attitude 
accuracy  such  that  multiple-site  measurements  (say  2  to  5  closely  spaced  sites)  might 
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achieve  the  improved  performance  suggested  by  the  results  shown  in  Figure  4.22  for 
the  case  of  no  uncertainty  in  position  and  attitude  for  a  2-site  data  set.  If  we  use  more 
sites,  can  we  relax  the  accuracy  standards? 
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4.3  Parameter  Extraction  Using  Multi-Point  Data  Sets  (Task  3) 


MMTargets  was  derived  from  ALLTEM-Double-Inversion,  which  was  designed  to  operate  on 
dynamically  acquired  data  from  the  USGS  ALLTEM  system.  MMTargets  simply  attempts  to 
invert  a  model  to  fit  the  dataset  it  is  directed  to  use.  This  can  be  a  single  static  snapshot,  multiple 
static  snapshots,  a  segment  from  a  single  survey  line  of  dynamically  acquired  data,  or  segments 
from  several  survey  lines  of  dynamically  acquired  data  (i.e.  a  patch  of  data).  The  USGS  used 
only  dynamic  data,  and  extracted  patches  of  data  that  were  small,  adjacent,  survey  line  segments. 

The  multi-point  inversion  capability  was  demonstrated  using  synthetic  data  in  Section  4.15 
(Multi-Site  Performance).  Here,  the  capability  is  demonstrated  by  combining  data  from  two 
survey  locations  at  Camp  Butner  that  were  presumed  to  be  two  distinct  targets  (see  Figure  4.23). 
After  combining  the  data  from  the  two  survey  points,  MMTargets  produces  the  polarizability 
curves  shown  in  Figure  4.24.  These  curves  clearly  suggest  that  the  target  is  rod-like;  and  this 
assertion  was  confirmed  when  a  37  mm  ordnance  was  recovered. 
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4.4  Multiple  Source  Detection  (Task  4) 

Data  from  targets  that  cannot  be  fit  using  simple  single-dipole  models  can  either  be  clutter  or 
multiple  targets.  Targets  with  a  small  volume  and  a  poor  fit  to  simple  dipole  models  are  likely 
clutter.  If  multiple  models  provide  similar  misfit  metrics,  this  may  indicate  that  a  complex  (non¬ 
dipole)  target  is  present  and  more  advance  analysis  algorithms  may  be  needed.  However,  targets 
with  large  volume  and  a  poor  fit  to  simple  dipole  models  may  be  clutter  or  multiple  objects  in 
close  proximity.  This  section  considers  the  latter  possibility.  Note  that  some  of  the  funds 
originally  allocated  for  this  task  were  redirected  towards  the  integration  of  the  previously 
discussed  routines  into  Oasis  montaj.  Subsequently,  the  scope  of  this  task  has  been  reduced  and 
testing  of  the  routines  in  this  section  has  not  moved  beyond  simulated  data. 

Characterization  and  classification  of  objects  in  close  proximity  is  a  difficult  task  that  is  often 
exacerbated  by  sparse  datasets.  A  single  cued  snapshot  does  not  have  sufficient  spatial  coverage 
and  density  for  effective  multi-target  characterization  and  classification.  This  is  illustrated  in  the 
figures  that  follow,  which  portray  the  behavior  modeled  fields  from  two  prolate  spheroid  targets 
in  close  proximity.  Table  4.2  contains  the  target  parameters  used  in  the  simulation.  In  these 
plots,  the  dots  mark  receiver  locations  and  the  background  colors  represent  field  values  at  the 
receiver  locations.  Figure  4.25  shows  synthetic  data  from  a  single  cued  snapshot,  and  Figure  4.26 
includes  additional  snapshots  that  provide  increased  data  density  to  help  resolve  the  complex 
high-frequency  nature  of  targets  in  close  proximity.  Note  that  some  of  the  peripheral  data  points 
in  Figure  4.26  still  have  a  relative  large  amount  of  energy. 
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Figure  4.25:  Plots  of  simulated  MetalMapper  data  for  two  oblate  spheroidal  targets  in  close  proximity  (targets  1  and  2 
from  Table  3.2)  The  dots  indicate  receiver  locations  from  a  single  cued  snapshot.  These  sparse  data  are  insufficient  to 
capture  the  high-frequency  character  of  the  anomalies. 
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Figure  4.26:  Plots  of  simulated  MetalMapper  data  for  two  oblate  spheroidal  targets  in  close  proximity.  Black  dots 
indicate  receiver  locations  from  four  cued  snapshots.  The  shaded  region  shows  the  footprint  of  the  MetalMapper 
cart.  These  data  are  sufficient  to  capture  the  high-frequency  character  of  the  anomalies. 


Table  4.2.  Target  parameters  used  to  generate  modeled  data  for  testing  frequency-domain  filters  and  multiple 
inversion  algorithms. 


Target 

o 

(S/m) 

hr 

X 

(m) 

y 

(m) 

z 

(m) 

Y axial 

(m) 

^ transverse 

(m) 

a 

(deg.  az.) 

P 

(deg.  inc.) 

i 

107 

150 

0.15 

-0.10 

-0.24 

0.451 

0.025 

10 

-23 

2 

107 

150 

-0.10 

0.17 

-0.20 

0.451 

0.025 

150 

73 

3 

107 

150 

0.6 

1.15 

-0.20 

0.451 

0.025 

0 

0,  90 

4 

107 

150 

0.4 

1.35 

-0.20 

0.451 

0.025 

0 

0,  90 

5 

107 

150 

0.1 

0.15 

-0.20 

0.451 

0.025 

30 

45 

6 

107 

150 

-0.1 

-0.15 

-0.20 

0.451 

0.025 

120 

-45 

4.4.1  Reduction-to-Pole  and  Downward  Continuation 

Song  et  al.  [7,  8]  have  shown  that  the  success  of  multi-target  inversions  is  critically  dependent  on 
high-quality  starting  models.  Song  attacked  this  problem  by  Monte  Carlo  generation  of  potential 
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starting  models,  and  using  those  models  whose  simulated  data  reasonably  fit  the  actual  data.  The 
present  work  attempts  to  reduce  this  time-consuming  process  by  selecting  starting  models  from 
filtered  datasets.  Two  filtering  methods  are  used:  reduction-to-pole  and  downward  continuation. 
The  goal  of  these  methods  is  to  reduce  the  interference  that  occurs  with  the  fields  generated  by 
objects  in  close  proximity,  and  to  provide  high-quality  starting  models  and  constraints  for  multi¬ 
target  inversion  routines.  The  benefit  of  the  reduction-to-pole  algorithm  is  that  it  centers  the 
anomaly  over  the  causative  body.  For  surveys  that  do  not  have  transmitter  coils  directly  over  the 
target,  this  method  will  center  the  anomaly  over  the  target.  The  downward  continuation  routine 
calculates  the  field  distribution  that  would  occur  if  the  survey  instrument  were  placed  at  a  level 
closer  to  the  target.  This  reduces  the  interference  between  fields  for  objects  in  close  proximity. 
The  downside  to  the  routines  is  that  they  can  generate  noise  and  artifacts.  Results  from  these 
routines  are  used  to  pick  high-quality  starting  models  for  the  multi-target  inversion  and  to 
provide  added  constraints  for  the  target’s  lateral  position. 

The  frequency-domain  operations  described  in  this  section  require  a  dataset  with  sufficient 
density  to  enable  Fourier  transforms.  Here  it  is  heuristically  demonstrated  that  two  adjacent 
passes  of  dynamic  MetalMapper  data  have  sufficient  data  density  to  employ  a  non-uniform 
discrete  Fourier  transform,  and  that  spatial  frequency-domain  techniques  can  be  employed. 
Figure  4.27  shows  simulated  data  plots  for  two  proximal  targets,  and  Figure  4.28  shows  the  same 
data  after  forward  and  reverse  non-uniform  discrete  Fourier  transforms.  Although  there  is  some 
distortion  of  the  original  data,  the  transformed  data  is  reasonably  close  to  the  original  data  and 
warrants  using  the  non-uniform  Fourier  transform  for  frequency-domain  processing.  Traditional 
Fourier  transform  techniques  require  data  on  a  uniform  grid.  For  non-uniform  survey  data,  one 
can  grid  the  data  onto  a  uniform  grid  before  invoking  a  traditional  Fourier  transform.  For  the 
problem  at  hand,  we  would  need  to  grid  both  the  calculated  transmitted  field  and  the  measured 
received  fields  for  all  coils.  This  would  markedly  increase  the  amount  of  data  to  be  processed, 
so  an  alternative  method  is  formulated  that  operates  directly  on  the  existing  data  with  a  non- 
uniform  Fourier  transform. 


Figure  4.27:  Left:  original  image  from  two  vertically  oriented  oblate  spheroid  targets  (target  3  and  4  from  Table  4.2). 
Right:  original  image  from  two  horizontally  oriented  oblate  spheroid  targets  (target  5  and  6  from  Table  4.2).  Colored 
dots  indicate  receiver  locations,  and  large  black  dots  (and  arrows)  indicate  true  target  positions.  These  maps  simulate 
data  acquired  along  two  adjacent  profiles  offset  by  a  distance  of  lm. 
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Figure  4.28:  Images  resulting  from  application  of  non-uniform  Fourier  transform  followed  by  uniform  inverse 
transform.  Left:  image  from  two  vertically  oriented  oblate  spheroid  targets  (target  3  and  4  from  Table  4.2).  Right: 
image  from  two  horizontally  oriented  oblate  spheroid  targets  (target  5  and  6  from  Table  4.2).  No  window  function  was 
used.  Colored  dots  indicate  data  points,  and  large  black  dots  (and  arrows)  indicate  true  target  position. 


Downward  continuation  is  a  process  that  attempts  to  calculate  the  fields  that  would  be  measured 
if  the  survey  plane  were  closer  to  the  target.  In  effect,  fields  are  calculated  that  would  have  been 
measured  if  the  survey  were  conducted  at  a  lower  elevation  that  was  closer  to  the  target.  Being 
closer  to  the  target  reduces  the  interference  between  fields  produced  by  objects  in  close 
proximity.  For  objects  oriented  at  an  oblique  angle,  the  measured  field  anomalies  can  be  skewed 
off-center  from  the  horizontal  target  location.  Downward  continuation  can  also  move  the 
anomaly  closer  to  the  true  horizontal  location  of  the  target. 

If  the  primary  field  does  not  have  vertical  polarization  when  energizing  a  target  (or  if  the  dipole 
magnetization  dips  at  an  oblique  angle),  then  the  measured  field  anomaly  will  be  skewed  off- 
center  from  the  horizontal  target  location.  Since  the  polarization  of  the  primary  field  at  the  target 
rotates  as  the  cart  moves,  many  polarizations  are  used  to  excite  the  target.  Reduction-to-pole  is  a 
process  that  calculates  the  anomaly  that  would  be  measured  if  the  primary  field  were  vertical. 
This  method  centers  the  measured  energy  over  the  target. 

This  section  provides  a  brief  mathematical  presentation  of  the  frequency-domain  routines,  and 
presents  results  that  illustrate  the  effectiveness  of  these  methods.  To  the  author’s  knowledge,  a 
reduction-to-pole  algorithm  suitable  for  EMI  data  has  not  been  previously  reported.  In  the 
frequency-domain  analysis  presented  below,  the  spatial  frequency  variables  are  obtained  from 
the  characteristic  equation,  with  Eigenvalues  A; ,  for  the  Laplace’s  equation.  See  Appendix  A  for 

more  background  on  the  notation  and  formulation  [9]: 

k2  =  Ar  +  Ar  +  Az  =  (ikx )  +  (/kv)  +  (kz)  =  (ikp )  +  k:  =  0  (4.10) 

where  kp  =  k2x  +  k2y ,  (x-y  in  horizontal  plane,  z  is  vertical)  which  leads  to 

K=\kp\-  (4.11) 

Beginning  with  the  field  of  a  dipole  target  in  the  spatial  domain: 
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the  fields  on  a  plane  above  the  target  at  the  origin  can  be  written  in  the  frequency-domain  so  the 
depth  (z  positive  up)  and  the  x-y  dependencies  are  separated: 
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A  dipole  source  at  an  arbitrary  location  is  given  by  phase  shifting  a  source  from  the  origin 

Barb  =  el(kxX+k’yy\kp\ zG0(k,z,m)  =  G(k,z,m,r') ,  (4.16) 

so  that  the  magnetization  of  an  arbitrary  collection  of  dipoles  is 

Banom  =  ^G(k,z,mn,rn).  (4.17) 

n 

If  the  induced  magnetization  by  the  primary  field  is  assumed  to  be  uniform,  than  the  single 
dipole  Green’s  function  can  be  factored  out: 


Banom=G0(k,z,m)^e 
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Reduction-to-pole  takes  the  form 
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and  downward  continuation  is  simply 
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(4.20) 


where  Az  is  the  downward  continuation  distance.  Finally,  both  of  these  operations  can  be 
combined 
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Conventional  reduction-to-pole  routines  usually  require  a  fixed  direction  for  the  primary  field. 
For  EMI  data,  the  orientation  of  the  primary  field  at  the  target  changes  with  cart  location  and 
transmitter  coil  orientation.  To  work  with  changing  primary  field  polarizations,  the  frequency- 
domain  phase  shift  in  equation  6.10  can  be  incorporated  into  the  Fourier  transform  used  to  move 
from  the  space-domain  to  the  frequency-domain.  The  discrete  Fourier  transform  projects  the 
space-domain  data  onto  frequency-domain  basis  functions,  and  each  sample  point  in  the  space- 
domain  has  a  frequency-domain  representation.  Since  the  primary  field  at  each  spatial  location 
is  different,  a  different  phase  shift  factor  multiplies  each  sample  point  in  the  frequency-domain 
before  it  is  summed: 


B(kx,ky)  =  ^ B(xt , yt )SRTPei('kxX‘+kyy‘) ,  (4.22) 

i 

where  Srtp  is  the  phase  shift  factor.  The 
changing  phase  shift  factor  is  applied  to  each 
data  point  in  the  space-domain  before 
calculating  the  Fourier  transform.  Changing  the 
order  of  operations  is  permissible  for  uniformly 
convergent  functions. 

In  practice,  both  reduction-to-pole  and 
downward  continuation  can  introduce  artifacts 
in  the  data,  and  can  produce  nonsensical  results 
if  misapplied.  Reduction-to-pole  has  well- 
known  difficulties  when  the  magnetization 
direction  approaches  horizontal.  Therefore, 
only  data  recorded  with  the  z-axis  transmitters 
and  receivers  are  used,  and  then  only  when  the 
primary  field  inclination  at  the  target  is  greater 
than  5  degrees.  The  difficulties  with  downward  continuation  arise  because  the  method  is 
essentially  a  high-pass  filter,  which  can  accentuate  high  frequency  noise  and  introduces  artifacts. 
Consider  that  the  spectrum  of  a  dipole  source  is  band-limited,  and  as  the  observation  point  moves 
farther  from  the  source  the  high  frequency  energy  is  reduced  (see  Figure  6.5).  Therefore,  the 
classical  dipole  spectrum  can  be  used  as  a  limiter  during  the  downward  continuation  process. 

The  reduction-to-pole  and  downward  continuation  methods  have  been  implemented  in  a  Python 
script.  Python  is  a  free  and  open-source  language  similar  to  Matlab  with  many  scientific 
libraries.  These  routines  were  applied  to  the  data  presented  in  Figure  4.27  Figure  4.28.  For 
calculating  the  primary  fields,  the  reduction-to-pole  routine  requires  a  target  position.  The 
centroid  of  the  anomaly  was  used  for  the  lateral  target  position,  and  0.25  m  was  used  for  the 
depth.  The  reduction-to-pole  results  are  shown  in  Figure  4.30  where  the  anomaly  peaks  have 
been  moved  very  close  to  the  actual  target  location  (compare  with  anomaly  peak  locations  in 
Figure  4.27  and  Figure  4.28).  Figure  4.31  shows  the  results  of  downward  continuation  for  10 
cm.  To  an  extent,  the  anomaly  energy  has  moved  closer  to  the  actual  target  locations,  but  some 
unwanted  artifacts  are  beginning  to  appear.  Figure  4.32  shows  the  combination  of  both 
reduction-to-pole  and  downward  continuation.  The  anomaly  peaks  are  very  close  to  the  actual 
target  positions,  but  there  are  many  unwanted  artifacts.  Low-pass  filtering  may  reduce  some  of 


Figure  4.29:  Frequency-domain  amplitude  spectra 
for  dipoles  at  various  depths.  The  amplitude  spectra 
are  independent  of  dipole  orientation. 
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these  artifacts,  but  noisy  data  will  cause  them  to  increase.  While  both  reduction-to-pole  and 
downward  continuation  each  help  to  provide  better  estimates  of  target  location,  the  best  results 
have  been  achieved  using  reduction-to-pole  with  no  downward  continuation.  These  preliminary 
results  suggest  that  more  accurate  starting  models  for  multi-target  inversion  can  be  obtained  with 
these  methods. 


Figure  4.30:  Images  resulting  from  the  reduction-to-pole  routine.  Left:  image  from  two  vertically  oriented  oblate 
spheroid  targets  (target  3  and  4  from  Table  6.1).  Right:  image  from  two  horizontally  oriented  oblate  spheroid 
targets  (target  5  and  6  from  Table  6.1).  Colored  dots  indicate  data  points,  and  large  black  dots  (and  arrows) 
indicate  true  target  positions.  Energy  peaks  are  now  very  close  to  actual  target  locations. 


Figure  4.31:  Images  resulting  from  application  downward  continuation  routines.  Left:  image  from  two  vertically 
oriented  oblate  spheroid  targets  (target  3  and  4  from  Table  3.2).  Right:  image  from  two  horizontally  oriented  oblate 
spheroid  targets  (target  5  and  6  from  Table  3.2).  Colored  dots  indicate  data  points,  and  large  black  dots  (and 
arrows)  indicate  true  target  positions. 


Gridding  the  data  before  reduction-to-pole  or  downward  continuation  may  reduce  artifacts,  but 
the  calculations  will  require  substantially  more  memory  and  longer  run  times.  Thus  far,  the 
routines  have  been  written  in  Python,  which  runs  slowly  compared  with  code  compiled  into 
native  machine  language  such  as  C  or  C++.  It  is  recommended  that  these  routines  be  rewritten  in 
a  faster  language. 
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image  from  two  horizontally  oriented  oblate  spheroid  targets  (target  5  and  6  from  Table  3.2).  Colored  dots 
indicate  data  points,  and  large  black  dots  (and  arrows)  indicate  true  target  positions.  Energy  peaks  are  very 
close  to  target  locations. 


4.4.2  Multi-Target  Inversion 

Inverting  for  multiple  targets  is  challenging  for  several  reasons.  First,  there  are  more  parameters 
to  solve  for.  Second,  all  of  the  parameters  have  a  non-linear  relationship  with  the  data  (when 
target  locations  are  unknown).  Once  the  locations  of  the  targets  are  known,  the  size  and 
orientation  can  be  determined  by  solving  a  linear  problem.  Nevertheless,  knowledge  of  the 
approximate  horizontal  target  position  can  significantly  improve  the  performance  of  multi-target 
inversion  routines. 

Several  groups  have  been  experimenting  with  multiple  target  solvers.  Song  et  al.  [7,8]  use  a 
Monte  Carlo  method  to  find  starting  models  with  a  reasonably  good  data  misfit,  and  then 
polished  the  parameter  estimates  using  conventional  decent  inversion  techniques.  Keiswetter 
[10]  uses  a  Monte  Carlo  technique  to  generate  an  initial  population  of  models,  then  selects  the 
best  combination  from  this  population,  and  finally  iterates  this  process  as  needed.  Shubitidze  et 
al.  [11]  characterizes  targets  using  ortho-normalized  volume  magnetic  sources.  Song  et  al.  [12] 
uses  the  MUSIC  algorithm  to  determine  multiple  target  positions,  and  then  applies  a  joint 
diagonalization  method  to  estimate  their  shape.  In  addition  to  MUSIC,  statistical  methods  such 
as  independent  component  analysis  and  blind  source  separation  can  be  applied  [13],  A 
combination  of  physical  constraints  such  as  downward  continuation  and  Green’s  function 
relationships  with  statistical  measures  such  as  independence  should  be  investigated. 

Because  this  project  was  multi-faceted,  the  limited  resources  for  a  particular  task  constrained 
most  developments  to  be  incremental  improvements  on  existing  routines.  Therefore,  for  the 
multi-target  inversion  we  examine  methods  that  did  not  require  coding  a  complex  routine  from 
scratch.  Rather,  we  attempted  to  extend  the  capabilities  of  the  existing  MMTargets  inversion 
routines.  This  involved  two  different  approaches.  The  first  is  a  simple  target  stripping  procedure 
and  the  second  is  the  interleaved  multi-target  solver. 
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The  simple  target  stripping  approach  was  investigated  by  creating  a  Python  script  to  invert  for 
two  targets,  one  at  a  time.  The  process  starts  by  constraining  the  lateral  position  of  the  targets 
based  on  the  reduction-to-pole  and  downward  continuation  methods  discussed  previously.  Then 
a  single  target  inversion  attempts  to  fit  the  data,  and  the  residual  dataset  is  calculated  by 
subtracting  the  modeled  data  for  the  best-fit  target.  Finally,  another  single  target  inversion  is  run 
to  fit  a  second  target  to  the  residual  data.  The  results  for  two  targets  in  close  proximity  at  the 
same  depth  were  poor.  When  trying  to  fit  one  target  to  a  multiple  target  dataset,  the  depths  are 
over  estimated  because  the  lateral  extent  of  multiple  targets  is  greater  than  individual  targets  and 


Table  4.3.  Single-target  inversion  results  (unconstrained  positions)  and  true  values  from  repeated  from  Table  6.1. 
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Points 

Target 
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(m) 

Y axial 

(m) 

^ transverse 

(m) 

a  (az.) 

(deg) 

P  (inc.) 

(deg) 

MSE 

(%) 

1  (true) 

0.15 

-0.10 

-0.24 

0.451 

0.025 

10 

-23 

2  (true) 

-0.10 

0.17 

-0.20 

0.451 

0.025 

150 

73 

i 

1+2 

-0.068 

0.171 

-0.276 

0.045 

0.044 

220 

-55.5 

3.7 

4 

1+2 

0.123 

0.345 

-0.275 

0.051 

0.037 

227 

65.5 

3.4 

5  (true) 

0.1 

0.15 

-0.20 

0.451 

0.025 

30 

45 

6  (true) 

-0.1 

-0.15 

-0.20 

0.451 

0.025 

120 

-45 

1 

5+6 

0.014 

0.010 

-0.301 

0.060 

0.041 

168 

-45 

3.3 

4 

5+6 

0.209 

0.210 

-0.285 

0.059 

0.038 

169 

-51 

3.4 

the  inversion  incorrectly  tries  to  fit  a  target  that  is  too  deep  to  better  match  the  data.  While  there 
may  be  some  cases  where  this  method  provides  reasonable  results,  the  performance  was 
unsatisfactory  for  cases  tested  in  this  report.  Therefore,  we  do  not  recommend  further  study  for 
this  method. 

The  interleaved  solver  takes  a  different  approach.  Because  the  cost  to  invert  a  dense  matrix 
scales  as  O (n2),  doubling  the  number  of  targets  increases  processing  time  by  a  factor  of  four  for 
each  iteration.  Rather  than  doubling  the  number  of  parameters  for  a  dual  target  system,  the 
MMTargets  multi-target  inversion  routine  considers  the  contribution  of  only  one  target  at  a  time. 
As  the  inversion  routine  iterates,  it  refines  the  parameter  estimates  for  one  target  at  a  time,  and 
repeatedly  cycles  through  the  entire  list  of  targets.  In  the  limit  of  iterative  parameter  updates 
shrinking  to  zero,  this  method  is  exactly  equivalent  to  a  simultaneous  multi-target  solver  using 
the  Gauss-Newton  method.  In  practice,  the  iterative  parameter  updates  are  small,  so  this  method 
reasonable  approximates  a  simultaneous  multi-target  solver.  With  this  method,  doubling  the 
number  of  targets  only  increases  the  solver  run  time  by  a  factor  of  two.  However,  results  show 
that  it  is  very  likely  that  iteratively  improving  the  parameter  estimates  for  a  given  set  of  starting 
models  will  reduce  the  data  misfit  to  a  local  minimum.  To  address  this  problem,  one  can  either 
specify  a  large  number  of  initial  models,  or  MMTargets  can  generate  these  models  randomly. 

To  evaluate  solver  performance,  simulated  data  from  target  pairs  from  Table  4.2  were  summed 
together.  Table  4.3  lists  solver  results  from  the  single  target  inversion.  For  single  target 


63 


Expanded  Processing  Techniques  for  EMI  Systems 


SERDP  MR- 1772 


inversions,  two  different  starting  models  were  used  with  essentially  no  constraints  on  parameter 
ranges.  For  the  multi-target  interleaved  solver,  two  initial  models  were  used  and  subsequent 
models  (100  in  total)  were  selected  randomly  (genetic  mutation  is  also  supported)  until  a  new 
model  with  a  mean  squared  error  (MSE)  less  than  90%  was  found.  The  proper  horizontal 
position  of  the  targets  was  given  in  the  initial  models  and  were  only  allowed  to  vary  by  +/-  1  cm. 
The  results  are  listed  in  Table  4.4  for  both  a  single  queued  snapshot  and  for  four  queued 
snapshots.  The  MSE  of  the  multi-target  solver  is  significantly  less  than  those  of  the  single  target 
solver.  The  estimated  depths  are  not  unreasonable,  but  vary  by  up  to  about  10  cm.  The  shape 
and  orientation  attributes  are  incorrect.  The  low  MSE  indicates  that  the  data  were  fitted 
reasonably  well  (real  data  with  misfits  of  2-3%  are  typical).  This  leads  to  the  conclusion  that  the 
multi-target  problem  is  ill-posed,  and  there  are  many  models  that  reasonably  fit  the  data. 

The  following  conclusions  and  recommendations  can  be  made  for  the  multiple  targets  problem 
based  on  the  test  cases  presented  earlier  in  this  Section. 

1.  Without  lateral  position  constraints,  the  interleaved  solver  does  not  converge  to  an 
acceptably  small  MSE.  Lateral  position  constraints  such  as  those  offered  by  the 
reduction-to-pole  routines  are  necessary  for  reasonable  multi-target  inversion 
performance. 

2.  The  multi-target  inversion  routines  and  frequency-domain  filters  have  only  been  tested 
with  modeled  data  with  no  noise.  A  much  more  thorough  evaluation  is  needed  using 
noisy  data  and  actual  survey  data. 

3.  Although  the  multi-source  inversion  problem  is  ill-posed,  the  routines  developed  for 
this  project  have  provided  reasonable  location  and  volume  estimates  for  the  cases 
tested. 

4.  Results  show  that  the  data  misfits  (MSE)  for  single  static  snapshots  are  similar  to  those 
from  multiple  queued  snapshots. 

5.  In  addition  to  MUSIC,  many  statistical  methods  such  as  independent  component 
analysis  and  blind  source  separation  may  be  leveraged  [13]  to  help  solve  the  multiple 
source  problem.  A  combination  of  physical  constraints  such  as  downward 
continuation  and  Green’s  function  relationships  with  statistical  measures  such  as 
independence  should  be  investigated. 

6.  At  each  Gauss-Newton  iteration,  the  interleaved  solver  attempts  to  perform  a  step-wise 
optimization  of  the  estimated  target  parameters  based  on  the  data  misfit.  For  the 
multi-target  problem,  local  minima  are  reached  very  frequently.  When  a  local 
minimum  is  reached,  a  new  set  of  parameters  is  created  via  Monte  Carlo.  Since  the 
Gauss-Newton  optimizations  oftentimes  does  not  make  substantial  progress  toward  a 
global  solution,  and  since  new  starting  models  are  found  by  Monte-Carlo,  the  routine 
is  essentially  a  Monte-Carlo  solver.  This  is  very  inefficient.  In  future  work,  we  plan 
to  develop  and  test  a  simultaneous  multiple-target  solver. 

7.  Choosing  new  starting  models  by  genetic  mutation  rather  than  Monte  Carlo  gave 
inferior  results. 
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8.  The  stripping  method  is  not  a  viable  method  for  multiple  target  inversion. 


Table  4.4.  Multi-target  inversion  results  (constrained  lateral  positions)  and  true  values  from  repeated  from  Table  6.1. 


Queued 

Points 

Target 

(m) 

y 

(m) 

Z 

(m) 

Y axial 

(m) 

^ transverse 

(m) 

a  (az.) 

(deg) 

P  (inc.) 

(deg) 

MSE 

(%) 

1  (true) 

0.15 

-0.10 

-0.24 

0.451 

0.025 

10 

-23 

2  (true) 

-0.10 

0.17 

-0.20 

0.451 

0.025 

150 

73 

i 

l:(l+2) 

0.140 

-0.090 

-0.227 

0.030 

0.030 
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60.7 

0.88 

i 

2:(l+2) 

-0.099 
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-0.184 

0.046 

0.024 

27.4 

-77.1 

0.88 

4 

l:(l+2) 

0.140 

-0.090 

-0.250 

0.051 

0.029 
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22.1 

0.59 

4 

2:(l+2) 

-0.110 

0.180 

-0.161 

0.041 

0.024 

156 

65.7 

0.59 

5  (true) 

0.1 

0.15 

-0.20 

0.451 

0.025 

30 

45 

6  (true) 

-0.1 

-0.15 

-0.20 

0.451 

0.025 

120 

-45 

1 

5:  (5+6) 

0.090 

0.150 

0.204 

0.047 

0.026 

358 

27.4 

1.3 

1 

6:  (5+6) 

-0.090 

-0.152 

-0.250 

0.039 

0.037 

239 

5.6 

1.3 

4 

5:  (5+6) 

0.090 

0.140 

-0.222 

0.047 

0.028 

202 

-35.6 

1.5 

4 

6:  (5+6) 

-0.090 

-0.140 

-0.250 
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0.033 

117 

-45.0 

1.5 
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4.5  Oasis  Montaj  Integration  (Task  5) 

The  task  of  integrating  some  of  the  routines  from  this  project  into  Geosoft’s  Oasis  montaj  was 
added  to  the  project  base  on  a  suggestion  from  the  program  office  during  an  IPR.  The  integrated 
routines  are  separated  into  two  groups.  This  first  group  contains  the  data  import  and  detection 
routines  for  dynamic  data.  These  routines  were  created  by  G  &  G  Geosciences,  Inc.,  and  are 
outlined  in  their  report.  The  second  group  contains  automated  MetalMapper  data  import  routines 
for  static  data,  and  automated  single  target  inversions  by  MMTargets.  This  group  of  routines  is 
available  to  Oasis  montaj  users  via  the  UX-Targets  menu.  There  are  four  menu  items  in  the  UX- 
Targets  menu: 

Convert  TEM  files  to  CSV  format 

Convert  TEM  noise  files  to  SU  format 


Import  ASCII  Target  Snapshot  Data 
Characterize  Targets 


The  remainder  of  this  section  presents  a  brief  overview  of  the  UX-Targets  routines.  A  UX- 
Targets  user  manual  is  available  with  more  detailed  information. 


The  first  three  menu 
items  are  used  to 
streamline 

importing  large 
MetalMapper 
datasets  into  Oasis. 
During  a  re¬ 
acquisition  survey,  a 
large  number  of 
TEM  files  will  be 
generated  -  one  or 
more  files  for  each 
target  in  the  re¬ 
acquisition  list. 
Importing  these  data 
is  a  three-step 
process.  First  the 
binary  TEM  files 
are  converted  to 


E^TEM  ->CSV 


Background  Filename  [”.teml 


QEI® 


[  Set  Bkgnd 


Coordinate  Conversion 
0  UTM 


Attitude  Correction 


Mag  Decl 
Platform  Ht 


mil 

I 


Output  File  Prefix 

Qbs  Data  Filefsl  f.teml 


Zone  |  10  Q 1 

GPS  Ant  Offset 

X  d°] 

Y  |  51 

Z  |  1.43 1  [ 

Output  Filename  (“CSV) 


□  Target  ID 


C: \w\demo  project\CalS  tatic00004. 
C: \w\demo  project\CalS  tatic00005. 
C:  \w\demo  project\CalS  taticOOOOG. 
C:  \w\demo  project\CalS  tatic00007. 
C:  \w\demo  project\CalS  tatic00008. 
C:  \w\demo  project\CalS  tatic00009. 
C:\w\demo  project\CalStatic00001 . 
C:\w\demo  projectSCalS  taticQ0Q02. 


[  Process  Selected  | 


|  Process  All  | 


Figure  4.33:  Dialog  for  converting 
TEM  files  to  CSV  files  to  import  into 
Oasis. 


TEM2Su  Convert 


SB® 


Bka/SU  Filename 


CAvAdemo  project\CalStatic00001  .tern 


Set  B kg/Proto  File 


groupBoxI 

□  Remove  Bkg  0  UTM  Cnvrt  UTM  Zone  10  " 


0  Single  Pt  Conversion 
Output  Filename 


Data  Pt  No 

□ 


C:\w\demo  project\CalS  taticOOOOl  .tern 


|  Cnvrt  Selected  | 

Cnvrt  All  | 

|  Load  SU  File  | 

Figure  4.34:  Dialog  for  converting  noise 
file  from  TEM  to  SU  format  for 
MMTargets . 


ASCII  so  they  can  be  imported  by  Oasis.  Invoking  the 
Convert  TEM  files  to  CSV  format  presents  the  dialog  shown  in 


Figure  4.33.  The  user  simply  drags  and  drops  all  of  the  files  to  be  converted  onto  the  dialog. 
Enter  the  other  information  requested  by  the  dialog  and  press  Process  All.  The  result  will 
be  a  group  CSV  files,  one  for  each  TEM  file.  Next,  MMTargets  requires  a  TEM  data  file  from  a 
target  free  location  so  that  it  can  estimate  the  EM  noise  level  in  the  data.  This  is  done  with  the 
Convert  TEM  noise  files  to  SU  format  menu  item.  Enter  the  files  names  in  the 
dialog  as  shown  in  Figure  4.34  and  press  Cnvrt  All.  In  the  final  step,  all  of  the  data  are 
imported  into  Oasis  using  the  Import  ASCII  Target  Snapshot  Data  menu  item  (see 
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Figure  4.35  and  Figure  4.36.)  Because  all  the  data  import  routines  operate  on  the  entire  list  of 
files  in  the  survey,  the  data  import  process  is  quite  efficient. 


Figure  4.35:  Selecting  CSV  files  for 
importing  into  Oasis . 


Figure  4.36:  Resulting  imported  data  after  selecting  CSV  files  for  importing 


into  Oasis. 


The  final  menu  item,  Characterize 
Targets,  starts  the  single  target  inversion 
process  (see  Figure  4.37).  The  inversion  routine 
is  run  for  all  selected  targets,  and  the  results  are 
placed  in  the  targets  database  (see  Figure  4.38). 

MMTargets  generates  a  rich  set  of  target 
properties  that  can  be  used  in  the  discrimination 
process,  which  are  listed  in  Table  4.5.  Also,  a 
summary  plot  as  shown  in  Figure  4.39  is 
generated  for  each  target  with  the  principal 
polarizability  curves  and  some  import  target 
parameter  estimates.  These  plots  can  be  quickly 
scanned  using  the  Microsoft  Windows  Picture 
Viewer  application.  Because  the  dialog  allows 
selection  of  all  of  the  queued  snapshots,  the  entire 
dataset  can  be  processed  with  very  little  human 
interaction.  MMTargets  has  a  large  number  of 
user-configurable  setting  that  control  how  the 
program  processes  and  inverts  each  dataset.  These  settings  are  configured  using  keyword 
parameters  defined  in  an  INI  file.  The  detailed  function  of  each  key  word  is  beyond  the  scope  of 
this  report.  For  more  details,  users  are  referred  to  the  MMTargets  user  manual. 


MMTargets  UXO  Characterize  Targets 


The  TEM  snapshot  database  should  be  open.  If  parameters  are  blank  they  will  be  prompted 
before  processing.  Otherwise  change  parameters  as  needed.  A  targets  database  will  be 
prompted  for  and  opened  during  processing. 


Source  database:  TEM-cued-data.gdb 
0  Current  Line 
0  Selected  Lines 
0  All  Lines 


Hardware  cal  file:  C : \w\GEOMETRICS-UXO-project  1  \conf ig\MMapper . ini  j  [  Select...  ] 

Noise  data  file:  C:\w\GEOMETRICS-UXO-projectl\Noise.SU  |  Select...  | 


Line:  9  of  9 


In . . . . . iniiiil 


Output  directory:  C:\...\Oasis  montaj\user\temp\MNITargetsChar-2011-10-21T17-29-32 


Plot  directory:  c:\w\geometrics-uxo-projectl\MMTargetsChar-2011-10-21Tl 7-29-32\ 


Figure  4.37:  Dialog  for  invoking  MMTargets  in 
single  target  inversion  mode. 
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Figure  4.38:  Example  of  a  targets  database  containing  results  from  MMTargets  inversion. 

Table  4.5:  MMTargets  single  target  properties  list. 


Target  Type 

Property 

Description 

Sphere 

x,  y,  z 

Target  location 

rl,r2,  r3 

Principle  axis  radii 

azimuth,  inclination 

Target  orientation 

MSE 

Mean  squared  error  between  the 
modeled  data  and  the  actual  data 

Oblate  spheroid 

x,  y,  z 

Target  location 

rl,r2,  r3 

Principle  axis  radii 

azimuth,  inclination 

Target  orientation 

MSE 

Mean  squared  error  between  the 
modeled  data  and  the  actual  data 

Prolate  spheroid 

x,  y,  z 

Target  location 

rl,r2,  r3 

Principle  axis  radii 

azimuth,  inclination 

Target  orientation 

MSE 

Mean  squared  error  between  the 
modeled  data  and  the  actual  data 

Ellipsoid 

x,  y,  z 

Target  location 

rl,r2,  r3 

Principle  axis  radii 

azimuth,  inclination 

Target  orientation 

MSE 

Mean  squared  error  between  the 
modeled  data  and  the  actual  data 

Best  overall 

x,  y,  z 

Target  location 

rl,r2,  r3 

Principle  axis  radii 

azimuth,  inclination 

Target  orientation 

MSE 

Mean  squared  error  between  the 
modeled  data  and  the  actual  data 
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3  target00003-MMTargets.png  -  Windows  Picture  and  Fax  Viewer 


0  IS® 


Figure  4.39:  Example  of  a  target  summary  plot. 
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5  Discussion  and  Recommendations 

This  project  has  addressed  a  number  of  important  issues  associated  with  UXO  remediation.  This 
Section  summarizes  the  contributions  of  this  project,  makes  conclusions  from  the  results 
obtained,  and  makes  recommendations  for  future  work. 

5.1  New  Contributions  and  Conclusions 

The  tasks  in  this  project  all  made  significant  contributions  to  UXO  research.  Before  considering 
the  contributions  of  each  task,  we  list  the  capabilities  of  a  program  named  MMTargets  that  has 
contributed  to  each  of  these  tasks.  MMTargets  has  been  configured  for  the  MetalMapper ,  and 
has  several  operating  modes  as  described  below. 

1 .  MMTargets  can  create  synthetic  datasets  to  simulate  both  static  and  dynamic  datasets  for 
a  wide  range  of  targets.  EM  noise,  positional  uncertainty,  and  orientation  uncertainty  can 
each  be  added  to  the  simulated  data.  This  is  very  useful  for  studies  pertaining  to 
optimizing  coil  geometries,  positioning  equipment,  and  survey  design. 

2.  MMTargets  can  conduct  shape-constrained  inversions. 

3.  MMTargets  provides  consistent  orientation  of  principal  polarizabilities  through  the  entire 
decay  curve. 

4.  MMTargets  can  perform  multi-target  inversions. 

5.  MMTargets  has  a  physics-based  detection  routine  that  helps  detect  targets  along  lines  of 
dynamically  acquired  data.  In  this  mode,  it  provides  initial  estimates  of  target  position 
and  size. 

6.  MMTargets  takes  coil  geometry  definitions  from  a  text  file,  so  it  can  easily  be  customized 
for  processing  data  acquired  by  other  time-domain  EMI  instruments. 


The  five  main  tasks  in  this  project  covered  a  lot  of  ground.  The  salient  contributions  of  each  task 
as  listed  below. 

1 .  Improved  Target  Detection  -  This  task  implements  a  simple  physics-based  indicator  in 
the  MetalMapper  acquisition  program  and  in  Oasis  montaj.  The  model  used  by  the 
indicator  is  a  spherical  target.  An  inversion  routine  using  this  simple  model  has  also  been 
added  to  MMTargets.  These  routines  will  help  users  find  the  location  of  targets  during 
acquisition  and  processing.  This  will  likely  reduce  the  number  of  statically  acquired 
targets  that  are  off-center. 

2.  Target  Parameter  Extraction  with  Multiple-Point  Static  or  Dynamic  Data  Sets  -  The 
capability  of  MMTargets  to  invert  multiple-point  static  datasets  for  single  and  multiple 
targets  has  been  demonstrated.  Both  modeled  data  and  actual  survey  data  were  used  in 
the  single  target  demonstration,  while  only  modeled  data  were  used  in  the  multi-target 
demonstration.  Essentially,  MMTargets  has  the  functionality  to  invert  any  MetalMapper 
dataset  acquired  by  any  means  for  single  or  multiple  targets  to  provide  characteristic 
information  to  a  classifier. 


70 


Expanded  Processing  Techniques  for  EMI  Systems 


SERDP  MR- 1772 


3.  Optimal  Array  Configuration(s)  -  The  optimal  array  task  addressed  the  issue  of  survey 
instrument  cost  and  complexity  versus  benefit  to  characterization  and  classification 
efforts.  Several  transmitter-receiver  permutations  were  investigated. 

4.  Multiple  Source  Detection  -  The  objective  here  is  to  identify  those  target  picks  associated 
with  complex  secondary  magnetic  fields  that  arise  either  from  the  presence  of  multiple 
targets  in  close  proximity,  targets  that  are  asymmetrical,  or  that  may  be  too  large 
dimensionally  to  be  approximated  with  a  point  dipole.  A  multi-target  inversion  algorithm 
has  been  developed  and  tested.  Downward  continuation  and  reduction-to-pole  algorithms 
have  also  been  developed  to  aid  interpretation  and  to  provide  useful  starting  models  for 
the  multi-target  inversion. 

5.  Oasis  montaj  Integration  -  This  task  was  not  part  of  the  original  work  scope,  but  was 
added  after  the  program  office  suggested  that  the  new  routines  should  be  integrated  into 
existing  Oasis  montaj  workflows.  TEM2CSV  and  TEM2SU  were  created  to  help  import 
the  TEM  files  recorded  by  the  MetalMapper  system  into  Oasis  montj.  Oasis  montaj  GXs 
have  been  written  to  do  the  following: 

a.  Streamline  the  import  of  MetalMapper  data, 

b.  Perform  constrained  inversions  using  several  basic  shapes, 

c.  Generate  inversion  summary  plots  for  each  target, 

d.  Import  inversion  results  into  Oasis  montaj  for  use  in  classification  routines. 


5.2  Recommendations  and  Future  Work 

This  project  has  expanded  the  software  support  for  current  state-of-the-practice  MetalMapper 
surveys.  The  detection  and  inversion  phases  have  been  enhanced.  These  routines  have  been 
demonstrated  with  modeled  data  and  with  a  few  cases  from  actual  surveys.  Nonetheless,  more 
testing  is  needed  demonstrate  that  these  routines  should  be  part  of  standard  workflows. 

5.2.1  Recommendation  for  MMTargets 

1.  The  primary  reason  for  the  broad  F-Test  distributions  is  noisy  data,  which  produces  a 
noise  objective  function.  While  the  noise  is  small  enough  to  permit  reasonable  inversion 
results,  the  noise  causes  a  number  of  small  local  minima  in  the  same  vicinity  of  the 
reported  solution.  The  local  minima  in  the  vicinity  of  the  solution  all  have  similar  MSE 
values,  and  some  will  likely  have  lower  MSE  values  than  that  of  the  reported  solution.  It 
is  anticipated  that  smoothing  the  objective  function  will  reduce  these  difficulties  and 
reduce  the  breadth  of  the  F-Test  distributions.  This  is  recommended  as  future  work. 

5.2.2  Recommendation  for  Detection 

1.  The  real-time  detection  routines  should  reduce  the  errors  and  erroneous  results  due  to 
static  data  collection  at  significant  lateral  offsets.  This  needs  to  be  demonstrated. 

2.  In  addition  to  the  state-of-the-practice  workflow,  we  have  demonstrated  the  capability  to 
invert  dynamic  data.  This  capability  needs  to  be  further  tested  and  integrated  into  the 
standard  Oasis  montage  workflows.  With  successful  dynamic  surveys,  it  may  be  possible 
to  significantly  reduce  the  cost  of  geophysical  detection  and  classification  efforts  by 
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significantly  reducing  or  eliminating  the  queued  re-acquisition  phase.  We  recommend 
further  development  and  testing  for  dynamic  surveys. 

5.2.3  Recommendations  for  Optimal  Arrays 

1.  Our  target  sets  were  generated  with  a  single  noise  level  consistent  with  noise 
measurements  at  several  demonstration  sites.  During  analysis  we  discovered  that  we  did 
not  account  properly  for  the  noise  reduction  due  to  stack  measurements  (VN).  As  a 
result,  the  noise  levels  were  perhaps  overly  high.  The  high  noise  does  not  affect  the 
relative  rankings  of  the  arrays  since  data  sets  for  each  array  were  generated  with  identical 
noise  statistics.  However,  we  would  like  to  investigate  different  noise  levels  and  how 
they  affect  the  overall  results. 

2.  We  would  like  to  extend  the  multiple  site  study  in  2  ways.  We  would  like  to 
experimentally  determine  the  required  specification  for  both  position  and  attitude 
accuracy  such  that  multiple-site  measurements  (say  2  to  5  closely  spaced  sites)  might 
achieve  the  improved  performance.  If  we  use  more  sites,  can  we  relax  the  standards  for 
platform  position  and  attitude  accuracy? 

5.2.4  Recommendations  for  Multiple  Targets 

1.  Multi-target  inversion  remains  difficult.  While  we  have  made  significant  inroads,  there 
are  still  testing  and  algorithmic  improvements  that  are  needed.  Reduction-to-pole 
reduces  the  search  space  by  significantly  constraining  the  lateral  position  of  the  target. 

2.  There  are  many  statistical  source  separation  methods  available  that  may  be  applicable  to 
the  multiple  source  problem.  In  addition  to  MUSIC,  statistical  methods  such  as 
independent  component  analysis  and  blind  source  separation  may  be  applicable  [9],  A 
combination  of  physical  constraints  such  as  downward  continuation  and  Green’s  function 
relationships  with  statistical  measures  such  as  independence  should  be  investigated. 

3.  The  interleaved  multi -target  solver  may  perform  better  if  the  parameters  for  all  targets  are 
simultaneously  optimized  at  each  Gauss-Newton  iteration.  This  should  be  tested. 

5.2.5  Recommendations  for  Oasis  montaj  Enhancements 

1 .  MMTargets  has  the  ability  to  model  static  and  dynamic  data  with  various  coil  geometries, 
EM  noise,  and  uncertainty  in  cart  location  and  orientation  provides  a  cost-effect  method 
for  testing  new  algorithms  and  proposed  instrument  arrays.  This  capability  should  be 
added  to  Oasis  montaj  to  allow  practitioners  to  design  surveys  and  test  various  scenarios. 

2.  The  reduction-to-pole  and  interleaved  multi-target  inversion  routines  should  be  integrated 
into  the  UXTargets  processing  menu  so  that  they  can  be  further  evaluated  using  actual 
survey  data. 

3.  In  order  to  invert  dynamic  data,  data  patches  outlined  by  polygons  must  be  generated. 
The  USGS  ALLTEM  data  processing  package  has  this  capability.  These  routines  should 
be  extended  to  operate  with  MMTargets  and  added  to  the  UXTargets  processing  menu. 

4.  The  USGS  ALLTEM  data  processing  package  has  a  set  of  statistical  classification 
routines  based  on  the  highly  regarded  open-source  ‘R’  statistical  software  package. 
These  routines  should  be  added  to  the  UXT argets  menu. 
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5.3  Future  Plans 

We  will  continue  to  seek  funding  sources  from  both  the  commercial  and  government  sectors  for 
testing  and  enhancing  the  new  processing  routines.  The  algorithms  we  have  developed  can 
easily  be  tailored  to  support  other  popular  EM  instruments  such  as  the  NanoTEM,  the  BUD,  the 
AOL,  the  TEMTADS,  and  the  EM63.  We  will  seek  demonstration  opportunities  using  these 
instruments  as  well. 
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7  Appendix  A 

We  describe  the  problem  formulation  for  Laplace’s  equation  using  the  separation  of  variables 
technique  below 


Starting  with  Laplace’s  equation:  X~B  = 


d2B  d2B  d2B 

a7+a7+a7 


=  0 


If  we  assume  that  B  is  separable,  that  is  it  can  be  written  as  a  product  of  three  functions  of  one 
variable  such  as  B=X(x)Y(y)Z(z),  then 


d2XYZ  d2XYZ  d2XYZ  d2Y  v„d2Z  1  d2X  1  d2Y  1  d2Z 

- ^— + - t —  + - r —  =  YZ — -  +  XZ — t  +  XY — -  = - -  + - -  + - - 

dx2  dy 2  dz  dx2  dy2  dz  X  dx2  Y  dy 2  Z  dz 


Now  the  equation  is  separated  and  each  term  is  a  function  of  one  variable.  The  only  way  this 
equation  can  be  satisfied  for  all  values  of  x,  y,  and  z,  is  if  each  term  has  a  constant  value  and  if 
the  sum  of  these  terms  is  zero. 


1  d2X 


1  d2Y 


1  d2Z  .  .  ,  . 

— -j  =  K + \ + K  =  0 


X  dxz  Y  dyz  Z  dz 

The  equation  with  the  X’s  is  the  characteristic  equation  for  the  Laplace  equation,  and  the  )Cs  are 
called  the  Eigenvalues.  If  we  assume  a  solution  of  the  Laplace  equation  of  the  form 


B  =  B(/v+a>y+rv'  =  X0ea*xY0ea?yZ0ea*z . 


then 


X0  ,  Y0  d2ea>y  ,  Z0  _  X0ay^  ,  Y0azye  ”  Z0a2e° 


X  dx2  Y  dy2  Z  dz2 


X 


■  =  a:+a,+  a2  =  0 


Now  at  least  one  term  (one  of  the  a’s)  in  the  above  equation  must  be  negative  and  at  least  one 
must  be  positive.  Although  there  are  other  possibilities,  we  are  free  to  choose  the  following, 
which  satisfies  this  sign  constraint  for  all  values  of  kx,  ky,  and  kz: 

ax  =  ikx\a  =  iky;az  =  k: .  This  leads  to  a  solution  of  the  form 


B  =  B0e 


i(kxx+kyy)+kzz 


=  X0eik‘%eik>yZ0ek *  =  X(x)Y(y)Z(z) 


which  casts  the  solution  in  terms  of  the  Fourier  transform  of  B  with  respect  to  x  and  y.  The 
wavenumber  variable  l<-  can  be  written  in  terms  of  the  wavenumber  variables  k,  and  kv: 


K=K+k2 
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